Workflow for Ridgway’s Hawk Integrated Population Model
Brian W. Rolek
2024-10-18
Supplemental materials for: Rolek, B.W., McClure, CJW, Dunn, L., Curti, M., … Ridgway’s Hawk IPM and PVA
Contact information: rolek.brian@peregrinefund.org
Metadata, data, and scripts used in analyses can be found at https://github.com/The-Peregrine-Fund/XXXXX.
The full workflow below is visible as a html website at: https://the-peregrine-fund.github.io/XXXXX/.
A permanent archive and DOI is available at: https://zenodo.org/doi/XXXXX
1. Code for IPM and PVA
1.1 Submodels without Integration
library('nimble')
library('parallel')
library ('coda')
load("/bsuscratch/brianrolek/riha_ipm/data.rdata")
#**********************
#* Parameter descriptions
#**********************
# Mark-resight-recovery data
# Observations (po) = y
# 1 seen first-year (age=0, just before 1st b-day)
# 2 seen nonbreeder
# 3 seen breeder
# 4 not seen
# States (ps)
# 1 alive first-year
# 2 alive nonbreeder
# 3 alive breeder
# 4 dead
# Groups
# 1 wild-born
# 2 translocated and hacked
###################################################
# PARAMETERS
# phiFY: survival probability first year
# phiA: survival probability nonbreeders
# phiB: survival probability breeders
# psiFYB: recruitment probability from first-year to breeder
# psiAB: recruitment probability from nonbreeders to breeder
# psiBA: recruitment probability from breeder to nonbreeders
# pFY: resight probability first-years
# pA: resight probability nonbreeders
# pB: resight probability breeders
# lmu.prod: mean productivity (males and females) per territory on the log scale
# sds: standard deviations for multivariate normal random effects for time and site
# R: correlation coefficients for multivariate normal random effects for time and site
# lambda: population growth rate (derived)
# extinct: binary indicator of extirpation at a site
# gamma: coefficient of effect from nest treatments
# betas: coefficient of effect from translocations
# deltas: coefficient of effect from survey effort
# mus: overall means for survival, recruitment, and detections
# r and rr: "r" parameter for negative binomial distribution
# also described as omega in manuscript
#**********************
#* Model code
#**********************
mycode <- nimbleCode(
{
###################################################
# Priors and constraints
###################################################
# Survival, recruitment, and detection can be correlated
for (k in 1:8){
betas[k] ~ dunif(-20, 20) # prior for coefficients
} # k
for (j in 1:8){
for (s in 1:nsite){
lmus[j,s] <- logit(mus[j,s])
mus[j,s] ~ dbeta(1,1) # prior for means
}} # m population #s sex #h hacked
# Temporal random effects and correlations among all sites, synchrony
for (j in 1:p){ sds[j] ~ dexp(1) }# prior for temporal variation estimated using the multivariate normal distribution
R[1:p,1:p] <- t(Ustar[1:p,1:p]) %*% Ustar[1:p,1:p] # calculate rhos, correlation coefficients
Ustar[1:p,1:p] ~ dlkj_corr_cholesky(eta=1.1, p=p) # Ustar is the Cholesky decomposition of the correlation matrix
U[1:p,1:p] <- uppertri_mult_diag(Ustar[1:p, 1:p], sds[1:p])
# Multivariate normal for temporal variance
for (t in 1:nyr){
eps[1:p,t] ~ dmnorm(mu.zeroes[1:p],
cholesky = U[1:p, 1:p], prec_param = 0)
}
# Temporal random effects and correlations between sites
for (jj in 1:p2){ sds2[jj] ~ dexp(1) }# prior for temporal variation estimated using the multivariate normal distribution
R2[1:p2,1:p2] <- t(Ustar2[1:p2,1:p2]) %*% Ustar2[1:p2,1:p2] # calculate rhos, correlation coefficients
Ustar2[1:p2,1:p2] ~ dlkj_corr_cholesky(eta=1.1, p=p2) # Ustar is the Cholesky decomposition of the correlation matrix
U2[1:p2,1:p2] <- uppertri_mult_diag(Ustar2[1:p2, 1:p2], sds2[1:p2])
# multivariate normal for temporal and site variance
for (t in 1:nyr){
for (s in 1:nsite){
eta[1:p2,s,t] ~ dmnorm(mu.zeroes2[1:p2],
cholesky = U2[1:p2, 1:p2], prec_param = 0)
} } # s t
###############################
# Likelihood for productivity
###############################
# Priors for productivity
for (s in 1:nsite){
lmu.f[s] ~ dnorm(0, sd=5)
} # s
gamma ~ dunif(-20, 20)
rr ~ dexp(0.05)
# Productivity likelihood
for (k in 1:nnest){
f[k] ~ dnegbin(ppp[k], rr)
ppp[k] <- rr/(rr+mu.f[k])
log(mu.f[k]) <- lmu.f[site.nest[k]] +
gamma*treat.nest[k] +
eps[9, year.nest[k] ] +
eta[9, site.nest[k], year.nest[k] ]
} # k
# Derive yearly brood size for population model
# Need to reorder because nimble doesn't
# handle nonconsecutive indices
# yrind.nest is a matrix of indices for each site
for (t in 1:nyr){
for (s in 1:nsite){
for (xxx in 1:nest.end[t,s]){
fecmat[t,s,xxx] <- mu.f[ yrind.nest[xxx,t,s] ]
} # xxx
mn.f[t,s] <- mean( fecmat[t,s,1:nest.end[t,s]] )
}} # s t
# GOF for number of fledglings
for (k in 1:nnest){
f.obs[k] <- f[k] # observed counts
f.exp[k] <- mu.f[k] # expected counts adult breeder
f.rep[k] ~ dnegbin(ppp[k], rr) # expected counts
f.dssm.obs[k] <- abs( ( f.obs[k] - f.exp[k] ) / (f.obs[k]+0.001) )
f.dssm.rep[k] <- abs( ( f.rep[k] - f.exp[k] ) / (f.rep[k]+0.001) )
} # k
f.dmape.obs <- sum(f.dssm.obs[1:nnest])
f.dmape.rep <- sum(f.dssm.rep[1:nnest])
f.tvm.obs <- sd(brood[1:nnest])^2/mean(brood[1:nnest])
f.tvm.rep <- sd(f.rep[1:nnest])^2/mean(f.rep[1:nnest])
################################
# Likelihood for counts
################################
# Omitted. See IPM or PVA.
################################
# Likelihood for survival
################################
# Calculate yearly averages for sites for integration
for (t in 1:nyr){
for (s in 1:nsite){
for (xxxx in 1:surv.end[t,s]){
# Need to reorder because nimble doesn't
# handle nonconsecutive indices
# yrind.surv is a matrix of indices for each site
phiFY2[ t, s, xxxx] <- phiFY[ yrind.surv[xxxx,t,s], t]
phiA2[ t, s, xxxx] <- phiA[ yrind.surv[xxxx,t,s], t]
phiB2[ t, s, xxxx] <- phiB[ yrind.surv[xxxx,t,s], t]
psiFYB2[ t, s, xxxx] <- psiFYB[ yrind.surv[xxxx,t,s], t]
psiAB2[ t, s, xxxx] <- psiAB[ yrind.surv[xxxx,t,s], t]
psiBA2[ t, s, xxxx] <- psiBA[ yrind.surv[xxxx,t,s], t]
pA2[ t, s, xxxx] <- pA[ yrind.surv[xxxx,t,s], t]
pB2[ t, s, xxxx] <- pB[ yrind.surv[xxxx,t,s], t]
} # xxxx
mn.phiFY[ t, s] <- mean( phiFY2[ t, s, 1:surv.end[t,s] ] )
mn.phiA[ t, s] <- mean( phiA2[ t, s, 1:surv.end[t,s] ] )
mn.phiB[ t, s] <- mean( phiB2[ t, s, 1:surv.end[t,s] ] )
mn.psiFYB[ t, s] <- mean( psiFYB2[ t, s, 1:surv.end[t,s] ] )
mn.psiAB[ t, s] <- mean( psiAB2[ t, s, 1:surv.end[t,s] ] )
mn.psiBA[ t, s] <- mean( psiBA2[ t, s, 1:surv.end[t,s] ] )
mn.pA[ t, s] <- mean( pA2[ t, s, 1:surv.end[t,s] ] )
mn.pB[ t, s] <- mean( pB2[ t, s, 1:surv.end[t,s] ] )
}}
for (i in 1:nind){
for (t in 1:nyr){
#Survival
logit(phiFY[i,t]) <- eta[1, site[i,t],t] + eps[1,t] +
lmus[1, site[i,t]] + betas[1]*hacked[i] # first year
logit(phiA[i,t]) <- eta[2, site[i,t],t] + eps[2,t] +
lmus[2, site[i,t]] + betas[2]*hacked[i] # nonbreeder
logit(phiB[i,t]) <- eta[3, site[i,t],t] + eps[3,t] +
lmus[3, site[i,t]] + betas[3]*hacked[i] # breeder
#Recruitment
logit(psiFYB[i,t]) <- eta[4, site[i,t],t] + eps[4,t] +
lmus[4, site[i,t]] + betas[4]*hacked[i] # first year to breeder
logit(psiAB[i,t]) <- eta[5, site[i,t],t] + eps[5,t] +
lmus[5, site[i,t]] + betas[5]*hacked[i] # nonbreeder to breeder
logit(psiBA[i,t]) <- #eta[6, site[i,t],t] + eps[6,t] +
lmus[6, site[i,t]] #+ betas[6]*hacked[i] # breeder to nonbreeder
#Re-encounter
logit(pA[i,t]) <- eta[7, site[i,t],t] + eps[7,t] +
lmus[7, site[i,t]] + betas[7]*hacked[i] # resight of nonbreeders
logit(pB[i,t]) <- eta[8, site[i,t],t] + eps[8,t] +
lmus[8, site[i,t]] + betas[8]*hacked[i] # resight of breeders
}#t
}#i
# Define state-transition and observation matrices
for (i in 1:nind){
# Define probabilities of state S(t+1) given S(t)
for (t in first[i]:(nyr-1)){
ps[1,i,t,1] <- 0
ps[1,i,t,2] <- phiFY[i,t] * (1-psiFYB[i,t])
ps[1,i,t,3] <- phiFY[i,t] * psiFYB[i,t]
ps[1,i,t,4] <- (1-phiFY[i,t])
ps[2,i,t,1] <- 0
ps[2,i,t,2] <- phiA[i,t] * (1-psiAB[i,t])
ps[2,i,t,3] <- phiA[i,t] * psiAB[i,t]
ps[2,i,t,4] <- (1-phiA[i,t])
ps[3,i,t,1] <- 0
ps[3,i,t,2] <- phiB[i,t] * psiBA[i,t]
ps[3,i,t,3] <- phiB[i,t] * (1-psiBA[i,t])
ps[3,i,t,4] <- (1-phiB[i,t])
ps[4,i,t,1] <- 0
ps[4,i,t,2] <- 0
ps[4,i,t,3] <- 0
ps[4,i,t,4] <- 1
# Define probabilities of O(t) given S(t)
po[1,i,t,1] <- 1
po[1,i,t,2] <- 0
po[1,i,t,3] <- 0
po[1,i,t,4] <- 0
po[2,i,t,1] <- 0
po[2,i,t,2] <- pA[i,t]
po[2,i,t,3] <- 0
po[2,i,t,4] <- (1-pA[i,t])
po[3,i,t,1] <- 0
po[3,i,t,2] <- 0
po[3,i,t,3] <- pB[i,t]
po[3,i,t,4] <- (1-pB[i,t])
po[4,i,t,1] <- 0
po[4,i,t,2] <- 0
po[4,i,t,3] <- 0
po[4,i,t,4] <- 1
} #t
} #i
# Likelihood
for (i in 1:nind){
# Define latent state at first capture
z[i,first[i]] <- y[i,first[i]]
for (t in (first[i]+1):nyr){
# State process: draw S(t) given S(t-1)
z[i,t] ~ dcat(ps[z[i,t-1], i, t-1, 1:4])
# Observation process: draw O(t) given S(t)
y[i,t] ~ dcat(po[z[i,t], i, t-1, 1:4])
} #t
} #i
} )
#**********************
#* Function to run model in NIMBLE
#**********************
run_ipm <- function(info, datl, constl, code){
library('nimble')
library('coda')
# helper function for multivariate normal
uppertri_mult_diag <- nimbleFunction(
run = function(mat = double(2), vec = double(1)) {
returnType(double(2))
p <- length(vec)
out <- matrix(nrow = p, ncol = p, init = FALSE)
for(i in 1:p){
out[ ,i] <- mat[ ,i] * vec[i]
}
return(out)
})
assign('uppertri_mult_diag', uppertri_mult_diag, envir = .GlobalEnv)
params <- c(# pop growth
#"lambda",
# fecundity
"lmu.f", "gamma", "rr", "mn.f",
# survival
"mus", "lmus", "betas",
# abundance
# "NB", "NF", "NFY", "N",
# "r",
# "N", "Ntot",
# error terms
"eps", "sds", "Ustar", "U", "R",
"eta", "sds2", "Ustar2", "U2", "R2",
# yearly summaries
'mn.phiFY', 'mn.phiA', 'mn.phiB',
'mn.psiFYB', 'mn.psiAB', 'mn.psiBA',
'mn.pA', 'mn.pB',
# goodness of fit
"f.dmape.obs", "f.dmape.rep",
"f.tvm.obs", "f.tvm.rep"
# "dmape.obs", "dmape.rep",
# "tvm.obs", "tvm.rep",
# "tturn.obs", "tturn.rep"
)
#n.chains=1; n.thin=200; n.iter=500000; n.burnin=300000
#n.chains=1; n.thin=50; n.iter=100000; n.burnin=50000
#n.chains=1; n.thin=10; n.iter=20000; n.burnin=10000
n.chains=1; n.thin=200; n.iter=500000; n.burnin=300000
mod <- nimbleModel(code,
constants = constl,
data = datl,
inits = info$inits,
buildDerivs = FALSE, # doesn't work when TRUE, no hope for HMC
calculate=T
)
cmod <- compileNimble(mod, showCompilerOutput = TRUE)
confhmc <- configureMCMC(mod)
confhmc$setMonitors(params)
hmc <- buildMCMC(confhmc)
chmc <- compileNimble(hmc, project = mod,
resetFunctions = TRUE,
showCompilerOutput = TRUE)
post <- runMCMC(chmc,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains,
thin = n.thin,
samplesAsCodaMCMC = T,
setSeed = info$seed)
return(post)
} # run_ipm function end
#*****************
#* Run chains in parallel
#*****************
this_cluster <- makeCluster(10)
post <- parLapply(cl = this_cluster,
X = par_info,
fun = run_ipm,
datl = datl,
constl = constl,
code = mycode)
stopCluster(this_cluster)
save(post, mycode,
file="/bsuscratch/brianrolek/riha_ipm/outputs/ipm_simp.rdata")
# save(post, mycode,
# file="C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\Ridgways IPM\\outputs\\ipm_statespace.rdata")1.2 Integrated Population Model
library('nimble')
library('parallel')
library ('coda')
load("/bsuscratch/brianrolek/riha_ipm/data.rdata")
#**********************
#* Parameter descriptions
#**********************
# Mark-resight-recovery data
# Observations (po) = y
# 1 seen first-year (age=0, just before 1st b-day)
# 2 seen nonbreeder
# 3 seen breeder
# 4 not seen
# States (ps)
# 1 alive first-year
# 2 alive nonbreeder
# 3 alive breeder
# 4 dead
# Groups
# 1 wild-born
# 2 translocated and hacked
###################################################
# PARAMETERS
# phiFY: survival probability first year
# phiA: survival probability nonbreeders
# phiB: survival probability breeders
# psiFYB: recruitment probability from first-year to breeder
# psiAB: recruitment probability from nonbreeders to breeder
# psiBA: recruitment probability from breeder to nonbreeders
# pFY: resight probability first-years
# pA: resight probability nonbreeders
# pB: resight probability breeders
# lmu.prod: mean productivity (males and females) per territory on the log scale
# sds: standard deviations for multivariate normal random effects for time and site
# R: correlation coefficients for multivariate normal random effects for time and site
# lambda: population growth rate (derived)
# extinct: binary indicator of extirpation at a site
# gamma: coefficient of effect from nest treatments
# betas: coefficient of effect from translocations
# deltas: coefficient of effect from survey effort
# mus: overall means for survival, recruitment, and detections
# r and rr: "r" parameter for negative binomial distribution
# also described as omega in manuscript
#**********************
#* Model code
#**********************
mycode <- nimbleCode(
{
###################################################
# Priors and constraints
###################################################
# survival, recruitment, and detection can be correlated
for (k in 1:8){
betas[k] ~ dunif(-20, 20) # prior for translocations coefficients
deltas[k] ~ dunif(-20, 20) # prior for survey effort coefficients
} # k
for (j in 1:8){
for (s in 1:nsite){
lmus[j,s] <- logit(mus[j,s])
mus[j,s] ~ dbeta(1,1) # prior for overall means
}} #
# Temporal random effects and correlations among all sites
for (j in 1:p){ sds[j] ~ dexp(1) }# prior for temporal variation
# estimated using the multivariate normal distribution
R[1:p,1:p] <- t(Ustar[1:p,1:p]) %*% Ustar[1:p,1:p] # calculate rhos, correlation coefficients
Ustar[1:p,1:p] ~ dlkj_corr_cholesky(eta=1.1, p=p) # Ustar is the Cholesky decomposition of the correlation matrix
U[1:p,1:p] <- uppertri_mult_diag(Ustar[1:p, 1:p], sds[1:p])
# multivariate normal for temporal variance
for (t in 1:nyr){ # survival params only have nyr-1, no problem to simulate from however
eps[1:p,t] ~ dmnorm(mu.zeroes[1:p],
cholesky = U[1:p, 1:p], prec_param = 0)
}
# Temporal random effects and correlations between sites
for (jj in 1:p2){ sds2[jj] ~ dexp(1) }# prior for temporal variation
# estimated using the multivariate normal distribution
R2[1:p2,1:p2] <- t(Ustar2[1:p2,1:p2]) %*% Ustar2[1:p2,1:p2] # calculate rhos, correlation coefficients
Ustar2[1:p2,1:p2] ~ dlkj_corr_cholesky(eta=1.1, p=p2) # Ustar is the Cholesky decomposition of the correlation matrix
U2[1:p2,1:p2] <- uppertri_mult_diag(Ustar2[1:p2, 1:p2], sds2[1:p2])
# multivariate normal for temporal variance
for (t in 1:nyr){ # survival params only have nyr-1, no problem to simulate from however
for (s in 1:nsite){
eta[1:p2,s,t] ~ dmnorm(mu.zeroes2[1:p2],
cholesky = U2[1:p2, 1:p2], prec_param = 0)
} } # s t
#######################
# Derived params
#######################
for (s in 1:nsite){
for (t in 1:(nyr-1)){
lambda[t, s] <- Ntot[t+1, s]/(Ntot[t, s])
loglambda[t, s] <- log(lambda[t, s])
}} #t
###############################
# Likelihood for productivity
###############################
# Priors
for (s in 1:nsite){
lmu.prod[s] ~ dnorm(0, sd=5)
} # s
gamma ~ dunif(-20, 20)
rr ~ dexp(0.05)
# Productivity likelihood
for (k in 1:npairsobs){
prod[k] ~ dnegbin(ppp[k], rr)
ppp[k] <- rr/(rr+mu.prod[k])
log(mu.prod[k]) <- lmu.prod[site.pair[k]] +
gamma*treat.pair[k] +
eps[9, year.pair[k] ] +
eta[9, site.pair[k], year.pair[k] ]
} # k
# Derive yearly productivity for population model
# need to reorder because nimble doesn't
# handle nonconsecutive indices
# yrind.pair is a matrix of indices for each site
for (t in 1:nyr){
for (s in 1:nsite){
for (xxx in 1:pair.end[t,s]){
prodmat[t,s,xxx] <- mu.prod[ yrind.pair[xxx,t,s] ]
} # xxx
mn.prod[t,s] <- mean( prodmat[t,s,1:pair.end[t,s]] )
}} # s t
# GOF for productivity
for (k in 1:npairsobs){
f.obs[k] <- prod[k] # observed counts
f.exp[k] <- mu.prod[k] # expected counts adult breeder
f.rep[k] ~ dnegbin(ppp[k], rr) # expected counts
f.dssm.obs[k] <- abs( ( f.obs[k] - f.exp[k] ) / (f.obs[k]+0.001) )
f.dssm.rep[k] <- abs( ( f.rep[k] - f.exp[k] ) / (f.rep[k]+0.001) )
} # k
f.dmape.obs <- sum(f.dssm.obs[1:npairsobs])
f.dmape.rep <- sum(f.dssm.rep[1:npairsobs])
f.tvm.obs <- sd(brood[1:npairsobs])^2/mean(brood[1:npairsobs])
f.tvm.rep <- sd(f.rep[1:npairsobs])^2/mean(f.rep[1:npairsobs])
################################
# Likelihood for counts
################################
# Abundance for year=1
for (v in 1:7){
for (s in 1:nsite){
# subtract one to allow dcat to include zero
N[v, 1, s] <- N2[v, 1, s] - 1
N2[v, 1, s] ~ dcat(pPrior[v, 1:s.end[v,s], s]) # Priors differ for FYs and adults
}} # s t
# Abundance for years > 1
for (t in 1:(nyr-1)){
for (s in 1:nsite){
# Number of wild born juvs
N[1, t+1, s] ~ dpois( (NFY[t, s]*mn.phiFY[t, s]*mn.psiFYB[t, s] + # first year breeders
NF[t, s]*mn.phiA[t, s]*mn.psiAB[t, s] + # nonbreeders to breeders
NB[t, s]*mn.phiB[t, s]*(1-mn.psiBA[t, s])) # breeders remaining
*mn.prod[t+1, s]/2 ) # end Poisson
# Abundance of nonbreeders
N[2, t+1, s] ~ dbin(mn.phiFY[t, s]*(1-mn.psiFYB[t, s]), NFY[t, s]) # Nestlings to nonbreeders
N[3, t+1, s] ~ dbin(mn.phiA[t, s]*(1-mn.psiAB[t, s]), NF[t, s]) # Nonbreeders to nonbreeders
N[4, t+1, s] ~ dbin(mn.phiB[t, s]*mn.psiBA[t, s], NB[t, s]) # Breeders to nonbreeders
# Abundance of breeders
N[5, t+1, s] ~ dbin(mn.phiFY[t, s]*mn.psiFYB[t, s], NFY[t, s]) # Nestlings to second year breeders
N[6, t+1, s] ~ dbin(mn.phiA[t, s]*mn.psiAB[t, s], NF[t, s]) # Nonbreeder to breeder
N[7, t+1, s] ~ dbin(mn.phiB[t, s]*(1-mn.psiBA[t, s]), NB[t, s]) # Breeder to breeder
}} # s t
# Count likelihoods, state-space model, and observation process
for (t in 1:nyr){
for (s in 1:nsite){
countsAdults[t, s] ~ dpois(lamAD[t, s]) # adult females
constraint_data[t, s] ~ dconstraint( (N[1, t, s] + hacked.counts[t, s]) >= 0 ) # constrain N1+hacked.counts to be >=0
NFY[t, s] <- N[1, t, s] + hacked.counts[t, s] # Transfers translocated first-year females
NF[t, s] <- sum(N[2:4, t, s]) # number of adult nonbreeder females
NB[t, s] <- sum(N[5:7, t, s]) # number of adult breeder females
NAD[t, s] <- NF[t, s] + NB[t, s] # number of adults
Ntot[t, s] <- sum(N[1:7, t, s]) # total number of females
log(lamAD[t, s]) <- log(NAD[t, s]) + deltas[1]*effort2[t, s] + deltas[2]*effort2[t, s]^2
log(lamFY[t, s]) <- log(N[1, t, s]) + deltas[3]*effort2[t, s] + deltas[4]*effort2[t, s]^2
}# s
# First-years at different sites have different distributions
# for better model fit
countsFY[t, 1] ~ dpois(lamFY[t, 1]) # doesn't have any zeroes so poisson
countsFY[t, 2] ~ dnegbin(pp[t], r) # first year females, includes translocated/hacked
pp[t] <- r/(r+(lamFY[t, 2] ))
} # t
r ~ dexp(0.05)
###################
# Assess GOF of the state-space models for counts
# Step 1: Compute statistic for observed data
# Step 2: Use discrepancy measure: mean absolute error
# Step 3: Use test statistic: number of turns
###################
for (t in 1:nyr){
c.repFY[t, 1] ~ dpois( lamFY[t, 1] )
c.repFY[t, 2] ~ dnegbin( pp[t], r )
for (s in 1:nsite){
c.expAD[t, s] <- lamAD[t, s] # expected counts adult breeder
c.expFY[t, s] <- lamFY[t, s]
c.obsAD[t, s] <- countsAdults[t, s]
c.obsFY[t, s] <- countsFY[t, s] # first year
c.repAD[t, s] ~ dpois( lamAD[t, s] ) # simulated counts
dssm.obsAD[t, s] <- abs( ( (c.obsAD[t, s]) - (c.expAD[t, s]) ) / (c.obsAD[t, s]+0.001) )
dssm.obsFY[t, s] <- abs( ( (c.obsFY[t, s]) - (c.expFY[t, s]) ) / (c.obsFY[t, s]+0.001) )
dssm.repAD[t, s] <- abs( ( (c.repAD[t, s]) - (c.expAD[t, s]) ) / (c.repAD[t, s]+0.001) )
dssm.repFY[t, s] <- abs( ( (c.repFY[t, s]) - (c.expFY[t, s]) ) / (c.repFY[t, s]+0.001) )
}} # t
dmape.obs[1] <- sum(dssm.obsAD[1:nyr, 1:nsite])
dmape.obs[2] <- sum(dssm.obsFY[1:nyr, 1:nsite])
dmape.rep[1] <- sum(dssm.repAD[1:nyr, 1:nsite])
dmape.rep[2] <- sum(dssm.repFY[1:nyr, 1:nsite])
tvm.obs[1] <- sd(c.obsB[1:nyr, 1:nsite])^2/mean(c.obsB[1:nyr, 1:nsite])
tvm.obs[2] <- sd(c.obsFY[1:nyr, 1:nsite])^2/mean(c.obsFY[1:nyr, 1:nsite])
tvm.rep[1] <- sd(c.repB[1:nyr, 1:nsite])^2/mean(c.repB[1:nyr, 1:nsite])
tvm.rep[2] <- sd(c.repFY[1:nyr, 1:nsite])^2/mean(c.repFY[1:nyr, 1:nsite])
# # Test statistic for number of turns
for (s in 1:nsite){
for (t in 1:(nyr-2)){
tt1.obsAD[t, s] <- step(c.obsAD[t+2, s] - c.obsAD[t+1, s])
tt2.obsAD[t, s] <- step(c.obsAD[t+1, s] - c.obsAD[t, s])
tt3.obsAD[t, s] <- equals(tt1.obsAD[t, s] + tt2.obsAD[t, s], 1)
tt1.obsFY[t, s] <- step(c.obsFY[t+2, s] - c.obsFY[t+1, s])
tt2.obsFY[t, s] <- step(c.obsFY[t+1, s] - c.obsFY[t, s])
tt3.obsFY[t, s] <- equals(tt1.obsFY[t, s] + tt2.obsFY[t, s], 1)
}} # t
tturn.obs[1] <- sum(tt3.obsAD[1:(nyr-2), 1:nsite])
tturn.obs[2] <- sum(tt3.obsFY[1:(nyr-2), 1:nsite])
for (s in 1:nsite){
for (t in 1:(nyr-2)){
tt1.repAD[t, s] <- step(c.repAD[t+2, s] - c.repAD[t+1, s])
tt2.repAD[t, s] <- step(c.repAD[t+1, s] - c.repAD[t, s])
tt3.repAD[t, s] <- equals(tt1.repAD[t, s] + tt2.repAD[t, s], 1)
tt1.repFY[t, s] <- step(c.repFY[t+2, s] - c.repFY[t+1, s])
tt2.repFY[t, s] <- step(c.repFY[t+1, s] - c.repFY[t, s])
tt3.repFY[t, s] <- equals(tt1.repFY[t, s] + tt2.repFY[t, s], 1)
}} # t
tturn.rep[1] <- sum(tt3.repAD[1:(nyr-2), 1:nsite])
tturn.rep[2] <- sum(tt3.repFY[1:(nyr-2), 1:nsite])
################################
# Likelihood for survival
################################
# Calculate averages for sites each year for integration
for (t in 1:nyr){
for (s in 1:nsite){
for (xxxx in 1:surv.end[t,s]){
# Reorder because nimble doesn't
# handle nonconsecutive indices
# yrind.surv is a matrix of indices for each site
phiFY2[ t, s, xxxx] <- phiFY[ yrind.surv[xxxx,t,s], t]
phiA2[ t, s, xxxx] <- phiA[ yrind.surv[xxxx,t,s], t]
phiB2[ t, s, xxxx] <- phiB[ yrind.surv[xxxx,t,s], t]
psiFYB2[ t, s, xxxx] <- psiFYB[ yrind.surv[xxxx,t,s], t]
psiAB2[ t, s, xxxx] <- psiAB[ yrind.surv[xxxx,t,s], t]
psiBA2[ t, s, xxxx] <- psiBA[ yrind.surv[xxxx,t,s], t]
pA2[ t, s, xxxx] <- pA[ yrind.surv[xxxx,t,s], t]
pB2[ t, s, xxxx] <- pB[ yrind.surv[xxxx,t,s], t]
} # xxxx
mn.phiFY[ t, s] <- mean( phiFY2[ t, s, 1:surv.end[t,s] ] )
mn.phiA[ t, s] <- mean( phiA2[ t, s, 1:surv.end[t,s] ] )
mn.phiB[ t, s] <- mean( phiB2[ t, s, 1:surv.end[t,s] ] )
mn.psiFYB[ t, s] <- mean( psiFYB2[ t, s, 1:surv.end[t,s] ] )
mn.psiAB[ t, s] <- mean( psiAB2[ t, s, 1:surv.end[t,s] ] )
mn.psiBA[ t, s] <- mean( psiBA2[ t, s, 1:surv.end[t,s] ] )
mn.pA[ t, s] <- mean( pA2[ t, s, 1:surv.end[t,s] ] )
mn.pB[ t, s] <- mean( pB2[ t, s, 1:surv.end[t,s] ] )
}}
for (i in 1:nind){
for (t in 1:nyr){
#Survival
logit(phiFY[i,t]) <- eta[1, site[i,t],t] + eps[1,t] +
lmus[1, site[i,t]] + betas[1]*hacked[i] # first year
logit(phiA[i,t]) <- eta[2, site[i,t],t] + eps[2,t] +
lmus[2, site[i,t]] + betas[2]*hacked[i] # nonbreeder
logit(phiB[i,t]) <- eta[3, site[i,t],t] + eps[3,t] +
lmus[3, site[i,t]] + betas[3]*hacked[i] # breeder
#Recruitment
logit(psiFYB[i,t]) <- eta[4, site[i,t],t] + eps[4,t] +
lmus[4, site[i,t]] + betas[4]*hacked[i] # first year to breeder
logit(psiAB[i,t]) <- eta[5, site[i,t],t] + eps[5,t] +
lmus[5, site[i,t]] + betas[5]*hacked[i] # nonbreeder to breeder
logit(psiBA[i,t]) <- #eta[6, site[i,t],t] + eps[6,t] +
lmus[6, site[i,t]] #+ betas[6]*hacked[i] # breeder to nonbreeder
#Re-encounter
logit(pA[i,t]) <- eta[7, site[i,t],t] + eps[7,t] +
lmus[7, site[i,t]] + betas[7]*hacked[i] +
deltas[5]*effort2[t, site[i,t]] + deltas[6]*effort2[t, site[i,t]]^2# resight of nonbreeders
logit(pB[i,t]) <- eta[8, site[i,t],t] + eps[8,t] +
lmus[8, site[i,t]] + betas[8]*hacked[i] +
deltas[7]*effort2[t, site[i,t]] + deltas[8]*effort2[t, site[i,t]]^2# resight of breeders
}#t
}#i
# Define state-transition and observation matrices
for (i in 1:nind){
# Define probabilities of state S(t+1) given S(t)
for (t in first[i]:(nyr-1)){
ps[1,i,t,1] <- 0
ps[1,i,t,2] <- phiFY[i,t] * (1-psiFYB[i,t])
ps[1,i,t,3] <- phiFY[i,t] * psiFYB[i,t]
ps[1,i,t,4] <- (1-phiFY[i,t])
ps[2,i,t,1] <- 0
ps[2,i,t,2] <- phiA[i,t] * (1-psiAB[i,t])
ps[2,i,t,3] <- phiA[i,t] * psiAB[i,t]
ps[2,i,t,4] <- (1-phiA[i,t])
ps[3,i,t,1] <- 0
ps[3,i,t,2] <- phiB[i,t] * psiBA[i,t]
ps[3,i,t,3] <- phiB[i,t] * (1-psiBA[i,t])
ps[3,i,t,4] <- (1-phiB[i,t])
ps[4,i,t,1] <- 0
ps[4,i,t,2] <- 0
ps[4,i,t,3] <- 0
ps[4,i,t,4] <- 1
# Define probabilities of O(t) given S(t)
po[1,i,t,1] <- 1
po[1,i,t,2] <- 0
po[1,i,t,3] <- 0
po[1,i,t,4] <- 0
po[2,i,t,1] <- 0
po[2,i,t,2] <- pA[i,t]
po[2,i,t,3] <- 0
po[2,i,t,4] <- (1-pA[i,t])
po[3,i,t,1] <- 0
po[3,i,t,2] <- 0
po[3,i,t,3] <- pB[i,t]
po[3,i,t,4] <- (1-pB[i,t])
po[4,i,t,1] <- 0
po[4,i,t,2] <- 0
po[4,i,t,3] <- 0
po[4,i,t,4] <- 1
} #t
} #i
# Likelihood
for (i in 1:nind){
# Define latent state at first capture
z[i,first[i]] <- y[i,first[i]]
for (t in (first[i]+1):nyr){
# State process: draw S(t) given S(t-1)
z[i,t] ~ dcat(ps[z[i,t-1], i, t-1, 1:4])
# Observation process: draw O(t) given S(t)
y[i,t] ~ dcat(po[z[i,t], i, t-1, 1:4])
} #t
} #i
} )
#**********************
#* Function to run model in NIMBLE
#**********************
run_ipm <- function(info, datl, constl, code){
library('nimble')
library('coda')
# helper function for multivariate normal
uppertri_mult_diag <- nimbleFunction(
run = function(mat = double(2), vec = double(1)) {
returnType(double(2))
p <- length(vec)
out <- matrix(nrow = p, ncol = p, init = FALSE)
for(i in 1:p){
out[ ,i] <- mat[ ,i] * vec[i]
}
return(out)
})
assign('uppertri_mult_diag', uppertri_mult_diag, envir = .GlobalEnv)
params <- c(
# pop growth
"lambda",
# fecundity
"lmu.prod", "gamma", "rr", "mn.prod",
# survival
"mus", "lmus", "betas", "deltas",
# abundance
"NB", "NF", "NFY", "N", "NAD",
"r",
"N", "Ntot",
# error terms
"eps", "sds", "Ustar", "U", "R",
"eta", "sds2", "Ustar2", "U2", "R2",
# yearly summaries
'mn.phiFY', 'mn.phiA', 'mn.phiB',
'mn.psiFYB', 'mn.psiAB', 'mn.psiBA',
'mn.pA', 'mn.pB',
# goodness of fit
"f.dmape.obs", "f.dmape.rep",
"f.tvm.obs", "f.tvm.rep",
"dmape.obs", "dmape.rep",
"tvm.obs", "tvm.rep",
"tturn.obs", "tturn.rep"
)
n.chains=1; n.thin=1; n.iter=500; n.burnin=100
#n.chains=1; n.thin=200; n.iter=600000; n.burnin=400000
mod <- nimbleModel(code,
constants = constl,
data = datl,
inits = info$inits,
buildDerivs = FALSE, # doesn't work when TRUE, no hope for HMC
calculate=T
)
cmod <- compileNimble(mod, showCompilerOutput = TRUE)
confhmc <- configureMCMC(mod)
confhmc$setMonitors(params)
hmc <- buildMCMC(confhmc)
chmc <- compileNimble(hmc, project = mod,
resetFunctions = TRUE,
showCompilerOutput = TRUE)
post <- runMCMC(chmc,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains,
thin = n.thin,
samplesAsCodaMCMC = T,
setSeed = info$seed)
return(post)
} # run_ipm function end
#*****************
#* Run chains in parallel
#*****************
this_cluster <- makeCluster(20)
post <- parLapply(cl = this_cluster,
X = par_info,
fun = run_ipm,
datl = datl,
constl = constl,
code = mycode)
stopCluster(this_cluster)
# save(post, mycode,
# file="/bsuscratch/brianrolek/riha_ipm/outputs/ipm_longrun.rdata")1.3 Population Viability Analysis
library('nimble')
library('parallel')
library ('coda')
load("/bsuscratch/brianrolek/riha_ipm/data.rdata")
source("/bsuscratch/brianrolek/riha_ipm/MCMCvis.R")
constl$K <- 100
cpus <- 10
#**********************
#* Extend "data" needed for PVA projections and scenarios
#**********************
constl$SC <- 12
hacked.counts <- array(0, dim=c(constl$SC, constl$nyr+constl$K, 3))
for (sc in 1:constl$SC){
hacked.counts[sc,1:constl$nyr,1:3] <- constl$hacked.counts
hacked.counts[sc,14:(constl$nyr+constl$K),1:3] <-
c(0, 0, 0, 1, 1, 1)[sc]*cbind(rep(-10, constl$K), rep(10, constl$K), rep(0, constl$K) )
}
constl$hacked.counts <- hacked.counts
datl$constraint_data <- rbind(datl$constraint_data, array(1, dim=c(constl$K,2)) )
#constl$treat.pair2 <- c(1, 0, 1, 1, 0, 1) # treated or not
constl$treat.pair2 <- c(0, 1, 1, 0, 1, 1) # treated or not
constl$hacked2 <- c(0, 0, 0, 1, 1, 1) # scenario- hacked or not
constl$effort2 <- rbind(constl$effort2, array(0, dim=c(constl$K,2), dimnames=list(2024:(2024+constl$K-1), c("LH", "PC"))))
# set number of treated nests for all scenarios
constl$num.treated <- c(0, 10, 100, 0, 10, 100, 0, 10, 100, 0, 10, 100) # 100s sub in for All but are over-ridden in model code
constl$treat.all <- c(0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1)
#**********************
#* Specify priors taken from posterior of IPM
#**********************
#* This helps ensure that all chains run
load("/bsuscratch/brianrolek/riha_ipm/outputs/ipm_longrun.rdata")
# Identify chains with NAs that failed to initialize
out <- lapply(post, as.mcmc)
NAlist <- c()
for (i in 1:length(out)){
NAlist[i] <- any (is.na(out[[i]][,1:286]) | out[[i]][,1:286]<0)
}
out <- out[!NAlist]
outp <- MCMCpstr(out, type="chains")
!NAlist
mn.lam <- apply(outp$lambda, c(2,3), mean)
lam <- apply(mn.lam, 1, mean)
Ni.func <- function (){
Ni <- outp$N[1:7,1:constl$nyr,1:2,
sample(seq(1, 5200, by=400), 1, replace = F)]
Ni.pva <- array(NA, dim=c(constl$SC, dim(Ni)+c(0,constl$K,0)))
for (sc in 1:constl$SC){
Ni.pva[sc, 1:7, 1:constl$nyr, 1:2] <- Ni
for (t in constl$nyr:(constl$nyr+constl$K)){
for (s in 1:2){
Ni.pva[sc, 1:7, t, s] <- Ni[1:7, 13, s] #round(Ni[1:7, 13, s]*lam[s]^(t-13))
}}} # t sc
return(Ni.pva)
} # function
u1 <- apply(outp$Ustar, c(1,2), mean)
u2 <- apply(outp$Ustar2, c(1,2), mean)
inits.func.pva <- function (){
list(
# fecundity
lmu.prod = apply(outp$lmu.prod, 1, mean),
gamma = mean(outp$gamma),
rr = mean(outp$rr),
# survival
z = par_info[[1]]$inits$z,
mus = apply(outp$mus, c(1,2), mean),
betas = apply(outp$betas, 1, mean),
deltas = apply(outp$deltas, 1, mean),
sds = apply(outp$sds, 1, mean),
Ustar = u1,
sds2 = apply(outp$sds2, 1, mean),
Ustar2 = u1,
# counts
countsAdults= matrix(c(374, 335, 305, 295, rep(NA, length(2015:2023)), rep(NA, length(2011:2023)) ), nrow=13),
r = mean(outp$r),
N = Ni.func()
)}
# Set seed then save for reproducibility
# then randomly draw for each chain
set.seed(1)
seeds <- sample(1:1000000, size=cpus, replace=FALSE)
par_info_pva <- list()
for (i in 1:cpus){
par_info_pva[[i]] <- list(seed=seeds[i], inits = inits.func.pva())
}
#**********************
#* Parameter descriptions
#**********************
# Mark-resight-recovery data
# Observations (po) = y
# 1 seen first-year (age=0, just before 1st b-day)
# 2 seen nonbreeder
# 3 seen breeder
# 4 not seen
# States (ps)
# 1 alive first-year
# 2 alive nonbreeder
# 3 alive breeder
# 4 dead
# Groups
# 1 wild-born
# 2 translocated and hacked
###################################################
# PARAMETERS
# phiFY: survival probability first year
# phiA: survival probability nonbreeders
# phiB: survival probability breeders
# psiFYB: recruitment probability from first-year to breeder
# psiAB: recruitment probability from nonbreeders to breeder
# psiBA: recruitment probability from breeder to nonbreeders
# pFY: resight probability first-years
# pA: resight probability nonbreeders
# pB: resight probability breeders
# lmu.prod: mean productivity (males and females) per territory on the log scale
# sds: standard deviations for multivariate normal random effects for time and site
# R: correlation coefficients for multivariate normal random effects for time and site
# lambda: population growth rate (derived)
# extinct: binary indicator of extirpation at a site
# gamma: coefficient of effect from nest treatments
# betas: coefficient of effect from translocations
# deltas: coefficient of effect from survey effort
# mus: overall means for survival, recruitment, and detections
# r and rr: "r" parameter for negative binomial distribution
# also described as omega in manuscript
#**********************
#* Model code
#**********************
mycode <- nimbleCode(
{
###################################################
# Priors and constraints
###################################################
# survival, recruitment, and detection can be correlated
for (k in 1:8){
betas[k] ~ dunif(-20, 20) # prior for coefficients
deltas[k] ~ dunif(-20, 20)
} # k
for (j in 1:8){
for (s in 1:nsite){
lmus[j,s] <- logit(mus[j,s])
mus[j,s] ~ dbeta(1,1) # prior for means
}} # m population #s sex #h hacked
# Temporal random effects and correlations among all sites
for (j in 1:p){ sds[j] ~ dexp(1) }# prior for temporal variation
# estimated using the multivariate normal distribution
R[1:p,1:p] <- t(Ustar[1:p,1:p]) %*% Ustar[1:p,1:p] # calculate rhos, correlation coefficients
Ustar[1:p,1:p] ~ dlkj_corr_cholesky(eta=1.1, p=p) # Ustar is the Cholesky decomposition of the correlation matrix
U[1:p,1:p] <- uppertri_mult_diag(Ustar[1:p, 1:p], sds[1:p])
# multivariate normal for temporal variance
for (t in 1:(nyr+K)){ # survival params only have nyr-1, no problem to simulate from however
eps[1:p,t] ~ dmnorm(mu.zeroes[1:p],
cholesky = U[1:p, 1:p], prec_param = 0)
}
# Temporal random effects and correlations between sites
for (jj in 1:p2){ sds2[jj] ~ dexp(1) }# prior for temporal variation
# estimated using the multivariate normal distribution
R2[1:p2,1:p2] <- t(Ustar2[1:p2,1:p2]) %*% Ustar2[1:p2,1:p2] # calculate rhos, correlation coefficients
Ustar2[1:p2,1:p2] ~ dlkj_corr_cholesky(eta=1.1, p=p2) # Ustar is the Cholesky decomposition of the correlation matrix
U2[1:p2,1:p2] <- uppertri_mult_diag(Ustar2[1:p2, 1:p2], sds2[1:p2])
# multivariate normal for temporal variance
for (t in 1:(nyr+K) ){ # survival params only have nyr-1, no problem to simulate from however
for (s in 1:nsite){
eta[1:p2,s,t] ~ dmnorm(mu.zeroes2[1:p2],
cholesky = U2[1:p2, 1:p2], prec_param = 0)
} } # s t
#######################
# Derived params
#######################
for(sc in 1:SC){
for (s in 1:nsite){
for (t in 1:(nyr-1+K)){
lambda[sc, t, s] <- Ntot[sc, t+1, s]/(Ntot[sc, t, s])
loglambda[sc, t, s] <- log(lambda[sc, t, s])
}} #t
for (s in 1:nsite){
for (t in 1:K){
extinct[sc, t, s] <- equals(Ntot[sc, nyr+t, s], 0)
extinctAD[sc, t, s] <- equals(NAD[sc, nyr+t, s], 0)
extinctB[sc, t, s] <- equals(NB[sc, nyr+t, s], 0)
}}
} # sc
###############################
# Likelihood for productivity
###############################
# Priors
for (s in 1:nsite){
lmu.prod[s] ~ dnorm(0, sd=5)
} # s
gamma ~ dunif(-20, 20)
rr ~ dexp(0.05)
# Productivity likelihood
for (k in 1:npairsobs){
prod[k] ~ dnegbin(ppp[k], rr)
ppp[k] <- rr/(rr+mu.prod[k])
log(mu.prod[k]) <- lmu.prod[site.pair[k]] +
gamma*treat.pair[k] +
eps[9, year.pair[k] ] +
eta[9, site.pair[k], year.pair[k] ]
} # k
# Derive yearly productivity for population model
# need to reorder because nimble doesn't
# handle nonconsecutive indices
# yrind.pair is a matrix of indices for each site
for (t in 1:nyr){
for (s in 1:nsite){
for (xxx in 1:pair.end[t,s]){
prodmat[t,s,xxx] <- mu.prod[ yrind.pair[xxx,t,s] ]
} # xxx
mn.prod[1,t,s] <- mean( prodmat[t,s,1:pair.end[t,s]] )
for(sc in 2:SC){
mn.prod[sc, t, s] <- mn.prod[1,t,s]
}}} # s t
# Future projections- fecundity
for(sc in 1:SC){
for (t in (nyr+1):(nyr+K) ){
for (s in 1:nsite){
# Calculate percent of nests treated
perc.treat[sc, t, s] <- treat.all[sc] + # If scenario =3 or 6, prop =1.0
step( NB[sc, t, s]-num.treated[sc] ) * (num.treated[sc] + 0.001)/(NB[sc, t, s]+0.001) * (1-treat.all[sc]) + # NB>=10, calculate proportion
step( NB[sc, t, s]-1 ) * 1-step( NB[sc, t, s]-num.treated[sc] ) * (1-treat.all[sc]) # NB>1 and NB<10 100%
log(mn.prod[sc, t, s]) <- lmu.prod[s] +
gamma*treat.pair2[sc]*perc.treat[sc, t, s] +
eps[9, t] +
eta[9, s, t]
}} } # s t sc
################################
# Likelihood for counts
################################
# Abundance for year=1
for (v in 1:7){
for (s in 1:nsite){
# Subtract one to allow dcat to include zero
N[1, v, 1, s] <- N2[1, v, 1, s] - 1
N2[1, v, 1, s] ~ dcat(pPrior[v, 1:s.end[v,s], s])
# Assign estimates for years with data to all scenarios
for (sc in 2:SC){
N[sc, v, 1, s] <- N[1, v, 1, s]
}
}} # s t
# Abundance for years > 1
for (t in 1:(nyr-1)){
for (s in 1:nsite){
# Number of wild born juvs
N[1, 1, t+1, s] ~ dpois( (NFY[1, t, s]*mn.phiFY[1, t, s]*mn.psiFYB[1, t, s] + # first year breeders
NF[1, t, s]*mn.phiA[1, t, s]*mn.psiAB[1, t, s] + # nonbreeders to breeders
NB[1, t, s]*mn.phiB[1, t, s]*(1-mn.psiBA[1, t, s])) # breeders remaining
*(mn.prod[1, t+1, s]/2) ) # end Poisson
# Abundance of nonbreeders
N[1, 2, t+1, s] ~ dbin(mn.phiFY[1, t, s]*(1-mn.psiFYB[1, t, s]), NFY[1, t, s]) # Nestlings to second year nonbreeders
N[1, 3, t+1, s] ~ dbin(mn.phiA[1, t, s]*(1-mn.psiAB[1, t, s]), NF[1, t, s]) # Nonbreeders to nonbreeders
N[1, 4, t+1, s] ~ dbin(mn.phiB[1, t, s]*mn.psiBA[1, t, s], NB[1, t, s]) # Breeders to nonbreeders
# Abundance of breeders
N[1, 5, t+1, s] ~ dbin(mn.phiFY[1, t, s]*mn.psiFYB[1, t, s], NFY[1, t, s]) # Nestlings to second year breeders
N[1, 6, t+1, s] ~ dbin(mn.phiA[1, t, s]*mn.psiAB[1, t, s], NF[1, t, s]) # Nonbreeder to breeder
N[1, 7, t+1, s] ~ dbin(mn.phiB[1, t, s]*(1-mn.psiBA[1, t, s]), NB[1, t, s]) # Breeder to breeder
# Assign estimates for years with data to
# all scenarios
for (v in 1:7){
for (sc in 2:SC){
N[sc, v, t+1, s] <- N[1, v, t+1, s]
} } # sc
}} # s t
# Future projections- population model
for (sc in 1:SC){
for (t in nyr:(nyr+K-1)){
for (s in 1:nsite){
# Number of wild born juvs
N[sc, 1, t+1, s] ~ dpois( (NFY[sc, t, s]*mn.phiFY[sc, t, s]*mn.psiFYB[sc, t, s] + # first year breeders
NF[sc, t, s]*mn.phiA[sc, t, s]*mn.psiAB[sc, t, s] + # nonbreeders to breeders
NB[sc, t, s]*mn.phiB[sc, t, s]*(1-mn.psiBA[sc, t, s])) # breeders remaining
*(mn.prod[sc, t+1, s]/2) ) # end Poisson
# Abundance of nonbreeders
## Second year nonbreeders
N[sc, 2, t+1, s] ~ dbin(mn.phiFY[sc, t, s]*(1-mn.psiFYB[sc, t, s]), NFY[sc, t, s]) # Nestlings to second year nonbreeders
## Adult nonbreeders
N[sc, 3, t+1, s] ~ dbin(mn.phiA[sc, t, s]*(1-mn.psiAB[sc, t, s]), NF[sc, t, s]) # Nonbreeders to nonbreeders
N[sc, 4, t+1, s] ~ dbin(mn.phiB[sc, t, s]*mn.psiBA[sc, t, s], NB[sc, t, s]) # Breeders to nonbreeders
# Abundance of breeders
## Second year breeders
N[sc, 5, t+1, s] ~ dbin(mn.phiFY[sc, t, s]*mn.psiFYB[sc, t, s], NFY[sc, t, s]) # Nestlings to second year breeders
## Adult breeders
N[sc, 6, t+1, s] ~ dbin(mn.phiA[sc, t, s]*mn.psiAB[sc, t, s], NF[sc, t, s]) # Nonbreeder to breeder
N[sc, 7, t+1, s] ~ dbin(mn.phiB[sc, t, s]*(1-mn.psiBA[sc, t, s]), NB[sc, t, s]) # Breeder to breeder
}} # s t
} # sc
# Count likelihoods, state-space model, and observation process
for (t in 1:nyr){
for (s in 1:nsite){
constraint_data[t, s] ~ dconstraint( (N[1, 1, t, s] + hacked.counts[1, t, s]) >= 0 ) # Transfers translocated first-year females
NFY[1, t, s] <- N[1, 1, t, s] + hacked.counts[1, t, s] # Transfers translocated first-year females
NF[1, t, s] <- sum(N[1, 2:4, t, s]) # number of adult nonbreeder females
NB[1, t, s] <- sum(N[1, 5:7, t, s]) # number of adult breeder females
NAD[1, t, s] <- NF[1, t, s] + NB[1, t, s] # number of adults
Ntot[1, t, s] <- sum(N[1, 1:7, t, s]) # total number of females
countsAdults[t, s] ~ dpois(lamAD[1, t, s]) # adult females
# constrain N1+hacked.counts to be >=0
log(lamAD[1, t, s]) <- log(NAD[1, t, s]) + deltas[1]*effort2[t, s] + deltas[2]*effort2[t, s]^2
log(lamFY[1, t, s]) <- log(N[1, 1, t, s]) + deltas[3]*effort2[t, s] + deltas[4]*effort2[t, s]^2
for (sc in 2:SC){
lamAD[sc, t, s] <- lamAD[1, t, s]
lamFY[sc, t, s] <- lamFY[1, t, s]
NFY[sc, t, s] <- NFY[1, t, s]
NF[sc, t, s] <- NF[1, t, s]
NB[sc, t, s] <- NB[1, t, s]
NAD[sc, t, s] <- NAD[1, t, s]
Ntot[sc, t, s] <- Ntot[1, t, s]
} # sc
}# s
# First-years at different sites have different distributions
# for better model fit
countsFY[t, 1] ~ dpois(lamFY[1, t, 1]) # doesn't have any zeroes so poisson
countsFY[t, 2] ~ dnegbin(pp[t], r) # first year females, includes translocated/hacked
pp[t] <- r/(r+(lamFY[1, t, 2] ))
} # t
r ~ dexp(0.05)
# Future projections of counts
for(sc in 1:SC){
for (t in (nyr+1):(nyr+K)){
for (s in 1:nsite){
# constrain N1+hacked.counts to be >=0, and allow it to shrink as the population shrinks
NFY[sc, t, s] <- step( (N[sc, 1, t, s] + hacked.counts[sc, t, s]) ) * (N[sc, 1, t, s] + hacked.counts[sc, t, s])
NF[sc, t, s] <- sum(N[sc, 2:4, t, s]) # number of adult nonbreeder females
NB[sc, t, s] <- sum(N[sc, 5:7, t, s]) # number of adult breeder females
NAD[sc, t, s] <- NF[sc, t, s] + NB[sc, t, s] # number of adults
Ntot[sc, t, s] <- sum(N[sc, 1:7, t, s]) # total number of females
}# s
} }# t sc
################################
# Likelihood for survival
################################
# Calculate yearly averages for integration
for (t in 1:(nyr-1)){
for (s in 1:nsite){
for (xxxx in 1:surv.end[t,s]){
# need to reorder because nimble doesn't
# handle nonconsecutive indices
# yrind.surv is a matrix of indices for each site
phiFY2[ t, s, xxxx] <- phiFY[ yrind.surv[xxxx,t,s], t]
phiA2[ t, s, xxxx] <- phiA[ yrind.surv[xxxx,t,s], t]
phiB2[ t, s, xxxx] <- phiB[ yrind.surv[xxxx,t,s], t]
psiFYB2[ t, s, xxxx] <- psiFYB[ yrind.surv[xxxx,t,s], t]
psiAB2[ t, s, xxxx] <- psiAB[ yrind.surv[xxxx,t,s], t]
psiBA2[ t, s, xxxx] <- psiBA[ yrind.surv[xxxx,t,s], t]
pA2[ t, s, xxxx] <- pA[ yrind.surv[xxxx,t,s], t]
pB2[ t, s, xxxx] <- pB[ yrind.surv[xxxx,t,s], t]
} # xxxx
mn.phiFY[1, t, s] <- mean( phiFY2[ t, s, 1:surv.end[t,s] ] )
mn.phiA[1, t, s] <- mean( phiA2[ t, s, 1:surv.end[t,s] ] )
mn.phiB[1, t, s] <- mean( phiB2[ t, s, 1:surv.end[t,s] ] )
mn.psiFYB[1, t, s] <- mean( psiFYB2[ t, s, 1:surv.end[t,s] ] )
mn.psiAB[1, t, s] <- mean( psiAB2[ t, s, 1:surv.end[t,s] ] )
mn.psiBA[1, t, s] <- mean( psiBA2[ t, s, 1:surv.end[t,s] ] )
mn.pA[1, t, s] <- mean( pA2[ t, s, 1:surv.end[t,s] ] )
mn.pB[1, t, s] <- mean( pB2[ t, s, 1:surv.end[t,s] ] )
for (sc in 2:SC){
mn.phiFY[sc, t, s] <- mn.phiFY[1, t, s]
mn.phiA[sc, t, s] <- mn.phiA[1, t, s]
mn.phiB[sc, t, s] <- mn.phiB[1, t, s]
mn.psiFYB[sc, t, s] <- mn.psiFYB[1, t, s]
mn.psiAB[sc, t, s] <- mn.psiAB[1, t, s]
mn.psiBA[sc, t, s] <- mn.psiBA[1, t, s]
mn.pA[sc, t, s] <- mn.pA[1, t, s]
mn.pB[sc, t, s] <- mn.pB[1, t, s]
}
}}
for (i in 1:nind){
for (t in 1:(nyr-1)){
#Survival
logit(phiFY[i,t]) <- eta[1, site[i,t],t] + eps[1,t] +
lmus[1, site[i,t]] + betas[1]*hacked[i] # first year
logit(phiA[i,t]) <- eta[2, site[i,t],t] + eps[2,t] +
lmus[2, site[i,t]] #+ betas[2]*hacked[i] # nonbreeder
logit(phiB[i,t]) <- eta[3, site[i,t],t] + eps[3,t] +
lmus[3, site[i,t]] + betas[3]*hacked[i] # breeder
#Recruitment
logit(psiFYB[i,t]) <- eta[4, site[i,t],t] + eps[4,t] +
lmus[4, site[i,t]] #+ betas[4]*hacked[i] # first year to breeder
logit(psiAB[i,t]) <- eta[5, site[i,t],t] + eps[5,t] +
lmus[5, site[i,t]] + betas[5]*hacked[i] # nonbreeder to breeder
logit(psiBA[i,t]) <- #eta[6, site[i,t],t] + eps[6,t] +
lmus[6, site[i,t]] #+ betas[6]*hacked[i] # breeder to nonbreeder
#Re-encounter
logit(pA[i,t]) <- eta[7, site[i,t],t] + eps[7,t] +
lmus[7, site[i,t]] + betas[7]*hacked[i] +
deltas[5]*effort2[t, site[i,t]] + deltas[6]*effort2[t, site[i,t]]^2 # resight of nonbreeders
logit(pB[i,t]) <- eta[8, site[i,t],t] + eps[8,t] +
lmus[8, site[i,t]] + #betas[8]*hacked[i] +
deltas[7]*effort2[t, site[i,t]] + deltas[8]*effort2[t, site[i,t]]^2 # resight of breeders
}#t
}#i
# Future projections of survival, recruitment, and detection
for(sc in 1:SC){
for (t in nyr:(nyr+K)){
for (s in 1:nsite){
# Calculate percent hacked
# assigns a one (1.0 or 100%) if abundance is less than 10
p.hackedFY[sc, t, s] <- step( NFY[sc, t, s]-10 ) * (10+0.001)/(NFY[sc, t, s]+0.001) + # NB>=10, calculate proportion
step( NFY[sc, t, s]-1 ) * 1-step( NFY[sc, t, s]-10 ) # NB>1 and NB<10 100%
p.hackedA[sc, t, s] <- step( NF[sc, t, s]-10 ) * (10+0.001)/(NF[sc, t, s]+0.001) + # NB>=10, calculate proportion
step( NF[sc, t, s]-1 ) * 1-step( NF[sc, t, s]-10 ) # NB>1 and NB<10 100%
p.hackedB[sc, t, s] <- step( NB[sc, t, s]-10 ) * (10+0.001)/(NB[sc, t, s]+0.001) + # NB>=10, calculate proportion
step( NB[sc, t, s]-1 ) * 1-step( NB[sc, t, s]-10 ) # NB>1 and NB<10 100%
# Hacking/translocation only affects birds
# moved to Punta Cana (site 2) hence- equals(s, 2)
#Survival
logit(mn.phiFY[sc, t, s]) <- eta[1, s, t] + eps[1, t] +
lmus[1, s] + betas[1]*hacked2[sc]*p.hackedFY[sc, t, s]*equals(s, 2) # first year
logit(mn.phiA[sc, t, s]) <- eta[2, s, t] + eps[2, t] +
lmus[2, s] #+ betas[2]*hacked2[sc]*p.hackedA[sc, t, s]*equals(s, 2) # nonbreeder
logit(mn.phiB[sc, t, s]) <- eta[3, s, t] + eps[3, t] +
lmus[3, s] + betas[3]*hacked2[sc]*p.hackedB[sc, t, s]*equals(s, 2) # breeder
#Recruitment
logit(mn.psiFYB[sc, t, s]) <- eta[4, s, t] + eps[4, t] +
lmus[4, s] #+ betas[4]*hacked2[sc]*p.hackedFY[sc, t, s]*equals(s, 2) # first year to breeder
logit(mn.psiAB[sc, t, s]) <- eta[5, s, t] + eps[5, t] +
lmus[5, s] + betas[5]*hacked2[sc]*p.hackedA[sc, t, s]*equals(s, 2) # nonbreeder to breeder
logit(mn.psiBA[sc, t, s]) <- #eta[6, site[i,t],t] + eps[6,t] +
lmus[6, s] #+ betas[6]*hacked[i] # breeder to nonbreeder
#Re-encounter
logit(mn.pA[sc, t, s]) <- eta[7, s, t] + eps[7, t] +
lmus[7, s] + betas[7]*hacked2[sc]*p.hackedA[sc, t, s]*equals(s, 2) # resight of nonbreeders
logit(mn.pB[sc, t, s]) <- eta[8, s, t] + eps[8, t] +
lmus[8, s] #+ betas[8]*hacked2[sc]*p.hackedB[sc, t, s]*equals(s, 2) # resight of breeders
} # s
} # t
} # sc
# Define state-transition and observation matrices
for (i in 1:nind){
# Define probabilities of state S(t+1) given S(t)
for (t in first[i]:(nyr-1)){
ps[1,i,t,1] <- 0
ps[1,i,t,2] <- phiFY[i,t] * (1-psiFYB[i,t])
ps[1,i,t,3] <- phiFY[i,t] * psiFYB[i,t]
ps[1,i,t,4] <- (1-phiFY[i,t])
ps[2,i,t,1] <- 0
ps[2,i,t,2] <- phiA[i,t] * (1-psiAB[i,t])
ps[2,i,t,3] <- phiA[i,t] * psiAB[i,t]
ps[2,i,t,4] <- (1-phiA[i,t])
ps[3,i,t,1] <- 0
ps[3,i,t,2] <- phiB[i,t] * psiBA[i,t]
ps[3,i,t,3] <- phiB[i,t] * (1-psiBA[i,t])
ps[3,i,t,4] <- (1-phiB[i,t])
ps[4,i,t,1] <- 0
ps[4,i,t,2] <- 0
ps[4,i,t,3] <- 0
ps[4,i,t,4] <- 1
# Define probabilities of O(t) given S(t)
po[1,i,t,1] <- 1
po[1,i,t,2] <- 0
po[1,i,t,3] <- 0
po[1,i,t,4] <- 0
po[2,i,t,1] <- 0
po[2,i,t,2] <- pA[i,t]
po[2,i,t,3] <- 0
po[2,i,t,4] <- (1-pA[i,t])
po[3,i,t,1] <- 0
po[3,i,t,2] <- 0
po[3,i,t,3] <- pB[i,t]
po[3,i,t,4] <- (1-pB[i,t])
po[4,i,t,1] <- 0
po[4,i,t,2] <- 0
po[4,i,t,3] <- 0
po[4,i,t,4] <- 1
} #t
} #i
# Likelihood
for (i in 1:nind){
# Define latent state at first capture
z[i,first[i]] <- y[i,first[i]]
for (t in (first[i]+1):nyr){
# State process: draw S(t) given S(t-1)
z[i,t] ~ dcat(ps[z[i,t-1], i, t-1, 1:4])
# Observation process: draw O(t) given S(t)
y[i,t] ~ dcat(po[z[i,t], i, t-1, 1:4])
} #t
} #i
} )
#**********************
#* Function to run model in NIMBLE
#**********************
run_pva <- function(info, datl, constl, code){
library('nimble')
library('coda')
# function for multivariate normal
uppertri_mult_diag <- nimbleFunction(
run = function(mat = double(2), vec = double(1)) {
returnType(double(2))
p <- length(vec)
out <- matrix(nrow = p, ncol = p, init = FALSE)
for(i in 1:p){
out[ ,i] <- mat[ ,i] * vec[i]
}
return(out)
})
assign('uppertri_mult_diag', uppertri_mult_diag, envir = .GlobalEnv)
params <- c(
# pop growth
"lambda",
# prod
"lmu.prod", "gamma", "rr", "mn.prod",
# survival
"mus", "lmus", "betas", "deltas",
# abundance
"NB", "NF", "NFY", "N", "NAD", "Ntot",
"r",
# error terms
"eps", "sds", "Ustar", "U", "R",
"eta", "sds2", "Ustar2", "U2", "R2",
# yearly summaries
'mn.phiFY', 'mn.phiA', 'mn.phiB',
'mn.psiFYB', 'mn.psiAB', 'mn.psiBA',
'mn.pA', 'mn.pB',
# pva
"extinct", "extinctAD", "extinctB"
)
#n.chains=1; n.thin=1; n.iter=500; n.burnin=100
n.chains=1; n.thin=200; n.iter=600000; n.burnin=400000
mod <- nimbleModel(code,
constants = constl,
data = datl,
inits = info$inits,
buildDerivs = FALSE,
calculate=T
)
cmod <- compileNimble(mod, showCompilerOutput = TRUE)
confhmc <- configureMCMC(mod)
confhmc$setMonitors(params)
hmc <- buildMCMC(confhmc)
chmc <- compileNimble(hmc, project = mod,
resetFunctions = TRUE,
showCompilerOutput = TRUE)
post <- runMCMC(chmc,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains,
thin = n.thin,
samplesAsCodaMCMC = T,
setSeed = info$seed)
return(post)
} # run_pva function end
#*****************
#* Run chains in parallel
#*****************
this_cluster <- makeCluster(cpus)
post <- parLapply(cl = this_cluster,
X = par_info_pva[1:cpus],
fun = run_pva,
datl = datl,
constl = constl,
code = mycode)
stopCluster(this_cluster)
save(post, mycode, seeds, cpus,
file="/bsuscratch/brianrolek/riha_ipm/outputs/pva.rdata")2. Figures and tables
2.1 Plot model estimates
library ('MCMCvis')
library ('coda')
library ('ggplot2')
library('reshape2')
library ("tidybayes")
library ('bayestestR')
library ("ggpubr")
library("viridis")
load("C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\Ridgways IPM\\outputs\\ipm_longrun.rdata")
load("data/data.rdata")
out <- lapply(post, as.mcmc)
outp <- MCMCpstr(out, type="chains")
#********************
#* Calculate summary stats
#********************
# Population growth rates
mn.lambda <- apply(outp$lambda, c(2,3), mean)
(md.lambda <- apply(mn.lambda, 1, median))## [1] 1.000060 1.304276(hdi95.lambda <- apply(mn.lambda, 1, HDInterval::hdi))## [,1] [,2]
## lower 0.9479604 1.209558
## upper 1.0348278 1.427075# Overall averages
mn.mus <- apply(outp$mus, c(1,3), mean)
(md.mus <- apply(mn.mus, 1, median))## [1] 0.33112222 0.93945299 0.88785788 0.11015826 0.40630345 0.07551559 0.20468205
## [8] 0.78345859(hdi.mus <- apply(mn.mus, 1, HDInterval::hdi))## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## lower 0.1721033 0.8122186 0.7967305 0.01096642 0.2414542 0.0327038 0.04307028
## upper 0.5341632 0.9994355 0.9647374 0.30829937 0.5460728 0.1214192 0.42815018
## [,8]
## lower 0.5467521
## upper 0.9503702# Compare between sites
df.sites <- data.frame(mus=1:8, pd=rep(NA,8), greater=NA)
for (i in 1:8){
first <- median(outp$mus[i, 1, ]) > median(outp$mus[i, 2, ])
df.sites$greater[i] <- ifelse(first, "LH", "PC")
if (first){
diff <- outp$mus[i, 1, ]-outp$mus[i, 2, ]
df.sites$pd[i] <- mean(diff>0) |> round(2)
} else{
diff <- outp$mus[i, 2, ]-outp$mus[i, 1, ]
df.sites$pd[i] <- mean(diff>0) |> round(2)
}
}
df.sites(md.mus <- apply(outp$mus, c(1,2), median)) |> round(2)## [,1] [,2]
## [1,] 0.39 0.26
## [2,] 0.97 0.93
## [3,] 0.91 0.87
## [4,] 0.12 0.07
## [5,] 0.11 0.70
## [6,] 0.12 0.03
## [7,] 0.06 0.34
## [8,] 0.71 0.89(hdi.mus <- apply(outp$mus, c(1,2), HDInterval::hdi)) |> round(2)## , , 1
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## lower 0.19 0.87 0.83 0.01 0.04 0.05 0.0 0.34
## upper 0.65 1.00 0.98 0.35 0.19 0.20 0.2 0.96
##
## , , 2
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## lower 0.07 0.68 0.74 0.0 0.41 0.00 0.04 0.6
## upper 0.52 1.00 0.97 0.4 0.95 0.08 0.74 1.0#*******************
#* Plot IPM results
#*******************
#Plot population size
Ntot <- outp$Ntot[1:13,,] # thin by 10 to speed up plotting
lN <- melt(Ntot)
colnames(lN) <- c("Time", "Site.ind", "Iter", "Abundance")
lN$Year <- lN$Time + 2010
lN$Site <- ifelse(lN$Site.ind==1, "Los Haitises", "Punta Cana")
med_hdis <- function(x, cm){
df <- data.frame(y=median(x),
ymin=HDInterval::hdi(x, credMass=cm)[1],
ymax=HDInterval::hdi(x, credMass=cm)[2] )
}
counts.tot <- datl$countsAdults+datl$countsFY+constl$hacked.counts[,-3]
df.counts <- melt(counts.tot)
colnames(df.counts) <- c("Year", "Site_ab", "Count")
df.counts$Site <- ifelse(df.counts$Site_ab=="LHNP", "Los Haitises", "Punta Cana")
p1 <- ggplot() + theme_minimal(base_size=14) +
geom_line(data=lN, aes(x=Year, y=Abundance, group=Iter),
color="gray30", linewidth=0.1, alpha=0.01 ) +
stat_summary(data=lN, aes(x=Year, y=Abundance),
geom="errorbar", width=0,
fun.data =med_hdis, fun.args=list(cm=0.95) ) +
stat_summary(data=lN, aes(x=Year, y=Abundance),
geom="errorbar", width=0, linewidth=1,
fun.data =med_hdis, fun.args=list(cm=0.80) ) +
stat_summary(data=lN, aes(x=Year, y=Abundance),
geom="line", linewidth=1,
fun.data =function(x){data.frame(y=median(x))} ) +
geom_line(data=df.counts, aes(x=Year, y=Count),
col="black", linewidth=1, linetype=2 ) +
facet_wrap("Site", scales="free_y") +
xlab("Year") + ylab("Number")
p1
# ggsave(filename="C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Abundance and counts.tiff",
# plot=p1,
# device="tiff",
# width=6, height=4, dpi=300)# Plot life stage structure
mdFY <- apply(outp$NFY[1:13,,], c(1,2), median)
mdB <- apply(outp$NB[1:13,,], c(1,2), median)
mdF <- apply(outp$NF[1:13,,], c(1,2), median)
lFY <- melt(mdFY)
lB <- melt(mdB)
lF <- melt(mdF)
lFY$Stage <- "First-year"
lB$Stage <- "Breeder"
lF$Stage <- "Nonbreeder"
ldat <- rbind(lFY, lB, lF)
colnames(ldat)[1:3] <- c("Year", "Sitenum", "Number")
ldat$Site <- ifelse(ldat$Sitenum==1, "Los Haitises", "Punta Cana")
ldat$year <- ldat$Year+2010
# Use median number of females in each stage
# to plot an approximate population structure
p3 <- ggplot(ldat, aes(fill=Stage, y=as.numeric(Number), x=year)) +
geom_bar(position="fill", stat="identity") +
ylab("Proportion of population") + xlab("Year") +
facet_wrap("Site") +
theme_minimal(base_size=14)+ theme(legend.position="top") +
scale_fill_viridis_d(option = "mako")
p4 <- ggplot(ldat, aes(fill=Stage, y=as.numeric(Number), x=year)) +
geom_bar(position="stack", stat="identity") +
ylab("Median number of females") + xlab("Year") +
facet_wrap("Site", scales = "free_y") +
theme_minimal(base_size=14 ) + theme(legend.position="top") +
scale_fill_viridis_d(option = "mako")
# plot with uncertainty
lFY <- melt(outp$NFY)
lB <- melt(outp$NB)
lF <- melt(outp$NF)
colnames(lFY) <- colnames(lB) <- colnames(lF) <- c("Time", "Site", "Iteration", "Abundance")
lFY$Year <- lFY$Time + 2010
lB$Year <- lB$Time + 2010
lF$Year <- lF$Time + 2010
lFY$Stage <- "First-year"
lB$Stage <- "Breeder"
lF$Stage <- "Nonbreeder"
lALL <- rbind(lFY, lB, lF)
lALL$Site <- ifelse(lALL$Site==1, "Los Haitises", "Punta Cana")
mdFY <- apply(outp$NFY[1:13,,], c(1,2), median)
mdB <- apply(outp$NB[1:13,,], c(1,2), median)
mdF <- apply(outp$NF[1:13,,], c(1,2), median)
hdiFY <- apply(outp$NFY[1:13,,], c(1,2), HDInterval::hdi)
hdiB <- apply(outp$NB[1:13,,], c(1,2), HDInterval::hdi)
hdiF <- apply(outp$NF[1:13,,], c(1,2), HDInterval::hdi)
medFY <- melt(mdFY)
medB <- melt(mdB)
medF <- melt(mdF)
hdisFY <- melt(hdiFY)
hdisB <- melt(hdiB)
hdisF <- melt(hdiF)
dfstages <- rbind(medFY, medB, medF)
dfhdis <- rbind(hdisFY, hdisB, hdisF)
dfstages$LHDI <- dfhdis[dfhdis$Var1=="lower",]
dfstages$UHDI <- dfhdis[dfhdis$Var1=="upper",]
cols <- mako(3)
med_hdis <- function(x, cm){
df <- data.frame(y=median(x),
ymin=HDInterval::hdi(x, credMass=cm)[1],
ymax=HDInterval::hdi(x, credMass=cm)[2] )
}
p5 <- ggplot() + theme_minimal(base_size=14) +
theme(legend.position="none") +
geom_line(data=lALL, aes(x=Year, y=Abundance,
group=Iteration, color=Stage,
alpha=Stage),
linewidth=0.1) +
stat_summary(data=lALL, aes(x=Year, y=Abundance),
geom="errorbar", width=0,
fun.data =med_hdis, fun.args=list(cm=0.95) ) +
stat_summary(data=lALL, aes(x=Year, y=Abundance),
geom="errorbar", width=0, linewidth=1,
fun.data =med_hdis, fun.args=list(cm=0.80) ) +
stat_summary(data=lALL, aes(x=Year, y=Abundance),
geom="line", linewidth=1,
fun.data =function(x){data.frame(y=median(x))} ) +
scale_color_manual( values = c("First-year"=cols[2], "Breeder"=cols[1], "Nonbreeder"=cols[3]) ) +
scale_alpha_manual( values = c("First-year"=0.015, "Breeder"=0.0075, "Nonbreeder"=0.15) ) +
facet_wrap(Stage~Site, scales="free_y",
shrink=TRUE, ncol=2) +
xlab("Year") + ylab("Number of females") +
theme(
strip.background = element_blank(),
strip.text.x = element_blank()
)
p35 <- ggarrange(p3, p5, ncol=1, nrow=2)
p35
# ggsave(filename="C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Stage structure_abundance.tiff",
# plot=p35,
# device="tiff",
# width=6, height=8, dpi=300)# plot survival, recruitment, and detection
mus.mat <- outp$mus
dimnames(mus.mat)[[1]] <- c("FY",
"A\nSurvival", "B",
"FY:B",
"A:B\nRecruitment",
"B:A",
"A", "B\nDetection")
dimnames(mus.mat)[[2]] <- c("Los Haitises", "Punta Cana")
lmus <- melt(mus.mat)
colnames(lmus) <- c("Parameter", "Site", "Iter", "value")
p6 <- ggplot(data= lmus, aes(x=value, y=Parameter )) +
stat_halfeye(point_interval=median_hdi, .width = c(.95, .80)) +
coord_flip() +
facet_wrap("Site") +
xlab("Probability") + ylab("Parameter") +
theme_bw(base_size=14) +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank())
p6
# ggsave(filename="C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Survival Recruitment Detection.tiff",
# plot=p5,
# device="tiff",
# width=8, height=4, dpi=300)
# Plot predicted survival, recruitment, and detection
# with effects from translocation and hacking
# Birds were only hacked from LHNP to PC here
# so we only predict values for PC
# pred.mus <- array(NA, dim=dim(outp$mus))
# for (m in 1:8){
# for (tr in 1:2){
# pred.mus[m,tr,] <- plogis( outp$lmus[m,2,] + outp$betas[m,]*c(0,1)[tr] )
# }}
#
# dimnames(pred.mus) <- list(param=c("FY", "A", "B", "FY:B", "A:B", "B:A", "A", "B"),
# trans=c("Not translocated", "Translocated"),
# iter=1:10000)
# lm <- melt(pred.mus)
#
# p5 <- ggplot(data= lm, aes(x=value, y=param )) +
# stat_halfeye(point_interval=median_hdi, .width = c(.95, .80)) +
# coord_flip() +
# facet_wrap("trans") +
# xlab("Probability") + ylab("Parameter") +
# theme_bw(base_size=14) +
# theme(panel.grid.major.x = element_blank(),
# panel.grid.minor.x = element_blank())
betas.pd <- apply(outp$betas, 1, pd)
ind <- which (betas.pd>0.975) # which are significant
pred.mus <- array(NA, dim=c(dim(outp$mus)))
for (m in ind){
for (tr in 1:2){
pred.mus[m,tr,] <- plogis( outp$lmus[m,2,] + outp$betas[m,]*c(0,1)[tr] )
}}
# treatment, mu, site, iter
mus.md <- apply(pred.mus, c(1,2), median)
mus.HDI80 <- apply(pred.mus, c(1,2), HDInterval::hdi, credMass=0.8)
mus.HDI95 <- apply(pred.mus, c(1,2), HDInterval::hdi)
# tiff(height=4, width=6, units="in", res=300,
# filename= "C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Translocation effects.tiff")
par(mar=c(3,4,4,1), mfrow=c(1,1))
for (tr in 1:2){
if (tr==1){
plot(1:4, c(0,0,1,1),
type="n",
pch=c(16, 17)[tr], cex=2,
ylim=c(0,1),
xlim=c(0.5, 4.5),
cex.axis=1, cex.lab=1.25,
ylab="Probability", xlab="",
main="",
xaxt="n", yaxt="n")
abline(h= seq(0,1,by=0.1), col="gray90")
abline(v= seq(0.5,6,by=1), col="gray90")
points(1:4+c(-0.1,0.1)[tr], mus.md[ind,tr],
pch=c(16, 17)[tr], cex=1.5)
}else{
points(1:4+c(-0.1,0.1)[tr], mus.md[ind,tr],
pch=c(16, 15)[tr], cex=1.5)
}}
axis(1, at=1:4, cex.axis=0.85,
labels=c("First-year\nSurvival",
"Breeder\nSurvival",
"Nonbreeder\n to Breeder\nRecruitment",
"Nonbreeder\nDetection"), las=1, padj=0.5)
axis(2, at=seq(0,1, by=0.1), cex.axis=1,
labels=c(0, NA,NA,NA,NA,0.5, NA,NA,NA,NA, 1))
for (m in 1:length(ind)){
for (tr in 1:2){
lines(rep(c(1:4)[m] + c(-0.1,0.1)[tr],2), mus.HDI95[1:2,ind[m],tr], lwd=2)
lines(rep(c(1:4)[m] + c(-0.1,0.1)[tr],2), mus.HDI80[1:2,ind[m],tr], lwd=4)
}}
legend(x=2.5,y=1.25, pch=c(16,15), pt.cex=1.5, cex=1, xpd=NA, horiz=T,
legend=c("Not translocated", "Translocated" ),
xjust=0.5)
#dev.off()# Plot fecundity in response to nest treatments
f <- exp(outp$lmu.prod)
f.md <- apply(f, 1, median)
f.HDI80 <- apply(f, 1, HDInterval::hdi, credMass=0.8)
f.HDI95 <- apply(f, 1, HDInterval::hdi)
# Calculate treatment effects on fecundity
f.pred <- array(NA, dim=dim(outp$lmu.prod))
for (s in 1:2){
f.pred[s,] <- exp(outp$lmu.prod[s,] + outp$gamma[,1])
} # s
f2.md <- apply(f.pred, 1, median)
f2.HDI80 <- apply(f.pred, 1, HDInterval::hdi, credMass=0.8)
f2.HDI95 <- apply(f.pred, 1, HDInterval::hdi)
# tiff(height=4, width=6, units="in", res=300,
# filename= "C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Fecundity_treatment.tiff")
par(mar=c(3,5,3,1), mfrow=c(1,1))
plot(c(1, 2)-0.1, f.md,
ylim=c(min(f.HDI95), max(f2.HDI95)),
xlim=c(0.5, 2.5),
ylab="Productivity", xlab="",
main="",
cex.axis=1.5, cex.lab=1.5,
xaxt="n", yaxt="n", type="n")
points(c(1, 2)-0.1, f.md, pch=16, cex=1.5)
abline(h=seq(0, 1.6, by=0.1), col="gray90")
abline(v=1.5, col="black")
axis(1, at=c(1,2), cex.axis=1,
labels=c("Los Haitises\nNational Park","Punta Cana"),
padj=0.5)
axis(2, at=seq(0, 1.4, by=0.1),
cex.axis=1,
labels=c(0, NA, NA, NA, 0.4, NA, NA, NA, 0.8, NA, NA, NA, 1.2, NA, NA))
lines(c(1,1)-0.1, f.HDI95[,1], lwd=2)
lines(c(2,2)-0.1, f.HDI95[,2], lwd=2)
lines(c(1,1)-0.1, f.HDI80[,1], lwd=4)
lines(c(2,2)-0.1, f.HDI80[,2], lwd=4)
points(c(1.1, 2.1), f2.md,
pch=18, cex=2.5)
lines(c(1,1) +0.1, f2.HDI95[,1], lwd=3)
lines(c(2,2) +0.1, f2.HDI95[,2], lwd=3)
lines(c(1,1) +0.1, f2.HDI80[,1], lwd=6)
lines(c(2,2) +0.1, f2.HDI80[,2], lwd=6)
legend(x=1.5,y=1.65, pch=c(16,18), pt.cex=c(1.5,2.5), cex=1,
legend=c("Untreated", "Treated" ),
horiz=TRUE, xjust=0.5, xpd=NA)
#dev.off()
# Calculate probability of direction
f.diff <- f[1,]-f[2,]
pd(f.diff)## Probability of Direction: 1.00f.md ## lmu.prod[1] lmu.prod[2]
## 0.3659702 0.2201483f.HDI95 ## lmu.prod[1] lmu.prod[2]
## lower 0.2877536 0.1581736
## upper 0.4512852 0.2860523f.diff2 <- f.pred[1,]-f.pred[2,]
pd(f.diff2)## Probability of Direction: 1.00f2.md ## [1] 1.1107340 0.6681587f2.HDI95## [,1] [,2]
## lower 0.8733436 0.4800632
## upper 1.3696686 0.8681799median(f.pred[1,]/f[1,])## [1] 3.03504median(f.pred[2,]/f[2,])## [1] 3.03504#*******************
# Correlations between demographic rates
#*******************
apply(outp$R, c(1,2), pd) ## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] 1.00000 0.51395 0.56615 0.61530 0.55105 0.50915 0.54825 0.50140 0.55795
## [2,] 0.51395 1.00000 0.50495 0.54665 0.54755 0.51475 0.58975 0.50035 0.51885
## [3,] 0.56615 0.50495 1.00000 0.58120 0.61845 0.51340 0.63790 0.51505 0.53490
## [4,] 0.61530 0.54665 0.58120 1.00000 0.60075 0.50615 0.71585 0.56350 0.68935
## [5,] 0.55105 0.54755 0.61845 0.60075 1.00000 0.50560 0.57150 0.54245 0.55390
## [6,] 0.50915 0.51475 0.51340 0.50615 0.50560 1.00000 0.51845 0.50560 0.50450
## [7,] 0.54825 0.58975 0.63790 0.71585 0.57150 0.51845 1.00000 0.52575 0.74405
## [8,] 0.50140 0.50035 0.51505 0.56350 0.54245 0.50560 0.52575 1.00000 0.53540
## [9,] 0.55795 0.51885 0.53490 0.68935 0.55390 0.50450 0.74405 0.53540 1.00000apply(outp$R2, c(1,2), pd) ## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] 1.00000 0.51055 0.51580 0.52700 0.53725 0.50115 0.50535 0.52995 0.52010
## [2,] 0.51055 1.00000 0.53865 0.50220 0.50155 0.50645 0.50730 0.56050 0.50395
## [3,] 0.51580 0.53865 1.00000 0.53130 0.54490 0.51285 0.55790 0.56290 0.52050
## [4,] 0.52700 0.50220 0.53130 1.00000 0.53680 0.50315 0.55635 0.65250 0.55840
## [5,] 0.53725 0.50155 0.54490 0.53680 1.00000 0.50265 0.52515 0.52370 0.55790
## [6,] 0.50115 0.50645 0.51285 0.50315 0.50265 1.00000 0.50840 0.51085 0.51740
## [7,] 0.50535 0.50730 0.55790 0.55635 0.52515 0.50840 1.00000 0.73010 0.52495
## [8,] 0.52995 0.56050 0.56290 0.65250 0.52370 0.51085 0.73010 1.00000 0.57520
## [9,] 0.52010 0.50395 0.52050 0.55840 0.55790 0.51740 0.52495 0.57520 1.00000median(outp$R2[7,8,]) ## [1] -0.1869805#*******************
# Correlations between demographics and population growth rates
#*******************
lam.m <- apply(outp$lambda[1:12,,], c(1,2), median)
lam.hdi <- apply(outp$lambda[1:12,,], c(1,2), HDInterval::hdi)
par(mfrow=c(1,2))
plot(2012:2023, lam.m[,1], type="b", pch=1,
ylab="Population growth rate", xlab="Year",
ylim=c(min(lam.hdi[,,1]), max(lam.hdi[,,1])),
main="Los Haitises")
abline(h=1, lty=2)
segments(x0=2012:2023, x1=2012:2023,
y0 = lam.hdi[1,,1], y1= lam.hdi[2,,1])
plot(2012:2023, lam.m[,2], type="b", pch=1, lty=1,
ylab="Population growth rate", xlab="Year",
ylim=c(min(lam.hdi[,,2]), max(lam.hdi[,,2])),
main="Punta Cana")
abline(h=1, lty=2)
segments(x0=2012:2023, x1=2012:2023,
y0 = lam.hdi[1,,2], y1= lam.hdi[2,,2])
# create a function to plot correlations
# between demographics and population growth rates
plot.cor <- function (lam.post, x.post, x.lab, ind.x=1:12, site=1, yaxt="b"){
# calculate the correlation coefficient
# over each iteration to propagate error
lambda <- lam.post[1:12, site,]
x <- x.post[ind.x, site, ]
df <- data.frame(ats1=c(0.8, 0.9, 1.0, 1.1, 1.2, 1.3),
ats2=c(1.0, 1.5, 2.0, 2.5, 3.0, 3.5),
labs1=c("0.8", NA, "1.0", NA, "1.2", NA),
labs2=c("1.0", NA, "2.0", NA, "3.0", NA)
)
if(site==1){ df <- df[, c(1,3)]} else{
df <- df[, c(2,4)]
}
lwd <- 1
cor.post <- array(NA, dim=dim(lambda)[-1])
for (i in 1:dim(lambda)[[2]]){
cor.post[i] <- cor(lambda[,i], x[,i])
}
cor.df <- data.frame(median= median(cor.post) |> round(digits=2),
ldi=HDInterval::hdi(cor.post)[1] |> round(digits=2),
hdi=HDInterval::hdi(cor.post)[2] |> round(digits=2),
pd = pd(cor.post) |> round(digits=2))
lam.m <- apply(lambda, 1, median)
lam.hdi <- apply(lambda, 1, HDInterval::hdi)
x.m <- apply(x, 1, median)
x.hdi <- apply(x, 1, HDInterval::hdi)
x.lims <- c(min(x.hdi[,]), max(x.hdi[,]))
y.lims <- c(min(lam.hdi[,]), max(lam.hdi[,]))
if(yaxt=="n"){
plot(x.m, lam.m[1:12],
xlim= x.lims,
ylim= y.lims,
type="n",
ylab="", xlab=x.lab, cex.lab=1.5,
main=c("", "")[site],
yaxt=yaxt, xaxt="n")
axis(2, at=df[,1], labels=c(rep(NA, 5)))
} else {
plot(x.m, lam.m[1:12],
xlim= x.lims,
ylim= y.lims,
type="n",
ylab="", xlab=x.lab, cex.lab=1.5,
main=c("", "")[site],
yaxt="n", xaxt="n")
axis(2, at=df[,1], labels=df[,2], las=1, cex.axis=1.5)
}
axis(1, cex.axis=1.5)
points(x.m, lam.m, pch=1)
segments(x0=x.hdi[1,], x1=x.hdi[2,],
y0 = lam.m, y1= lam.m, lwd=lwd)
segments(x0=x.m, x1=x.m,
y0 = lam.hdi[1,], y1= lam.hdi[2,], lwd=lwd)
xs <- c(x.lims[1], x.lims[1]+(x.lims[2]-x.lims[1]*1.5))
ys <- c( (y.lims[2]-y.lims[1])+y.lims[1], (y.lims[2]-y.lims[1])+y.lims[1],
(y.lims[2]-y.lims[1])*0.6+y.lims[1], (y.lims[2]-y.lims[1])*0.6+y.lims[1])
polygon(x=c(xs, rev(xs)), y=ys, col=alpha("white", 0.7), border=NA )
text(x = x.lims[1], y = (y.lims[2]-y.lims[1])*0.9+y.lims[1],
paste("r = ", cor.df$median, " (",cor.df$ldi,", ", cor.df$hdi, ")", sep=""),
pos = 4, font = 3, cex = 1.5)
text(x = x.lims[1], y = (y.lims[2]-y.lims[1])*0.7+y.lims[1], paste("P(r>0) = ", cor.df$pd, sep=""),
pos = 4, font = 3, cex = 1.5)
}
# Plor correlations between population grrowth rates
# and demographics.
# "r" is a correlation coefficient and represents
# the magnitude of the correlation
# P(r>0) is the probability of direction (similar to p-values)
# that is, the probability that an effect exists
# tiff(height=8, width=8, units="in", res=300,
# filename= "C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//PopGrowthand Demographics1.tiff")
# par(mfrow=c(4,3), mar=c(4,1,1,1), oma=c(0,5,3,0))
par(mfrow=c(2,3), mar=c(4,1,1,1), oma=c(0,5,3,0))
plot.cor(outp$lambda, outp$mn.prod, x.lab="Fecundity", ind.x=2:13)
plot.cor(outp$lambda, outp$mn.phiFY, x.lab="First-year Survival", yaxt="n")
plot.cor(outp$lambda, outp$mn.phiA, x.lab="Nonbreeder Survival", yaxt="n")
plot.cor(outp$lambda, outp$mn.phiB, x.lab="Breeder Survival")
plot.cor(outp$lambda, outp$mn.psiFYB, x.lab="First-year to Breeder", yaxt="n")
plot.cor(outp$lambda, outp$mn.psiAB, x.lab="Nonbreeder to Breeder", yaxt="n")
mtext("Population growth rate", side=2, outer=TRUE, line=2.5, cex=2)
mtext("Los Haitises", side=3, outer=TRUE, cex=2)
#dev.off()
# tiff(height=4, width=8, units="in", res=300,
# filename= "C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//PopGrowthand Demographics2.tiff")
par(mfrow=c(2,3), mar=c(4,1,1,1), oma=c(0,5,3,0))
plot.cor(outp$lambda, outp$mn.prod, x.lab="Fecundity", ind.x=2:13, site=2)
plot.cor(outp$lambda, outp$mn.phiFY, x.lab="First-year Survival", yaxt="n", site=2)
plot.cor(outp$lambda, outp$mn.phiA, x.lab="Nonbreeder Survival", yaxt="n", site=2)
plot.cor(outp$lambda, outp$mn.phiB, x.lab="Breeder Survival", site=2)
plot.cor(outp$lambda, outp$mn.psiFYB, x.lab="First-year to Breeder", yaxt="n", site=2)
plot.cor(outp$lambda, outp$mn.psiAB, x.lab="Nonbreeder to Breeder", yaxt="n", site=2)
mtext("Population growth rate", side=2, outer=TRUE, line=2.5, cex=2)
mtext("Punta Cana", side=3, outer=TRUE, cex=2)
#dev.off()#*******************
#* Plot PVA results
#*******************
# Abundance of females at Los Haitises
load("C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\Ridgways IPM\\outputs\\pva.rdata")
out <- lapply(post, as.mcmc)
outp <- MCMCpstr(out, type="chains")
library (viridis)
cols <- viridis(3)
pe <- apply(outp$extinct, c(1,2,3), mean)
# tiff(height=3, width=8, units="in", res=300,
# filename= "C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Quasiextinction.tiff")
par(mar=c(2, 2, 1, 1), oma=c(3,3,0,0))
layout(matrix(c(1,1,2,2,3), 1, 5, byrow = TRUE))
plot(2024:2123, pe[1, ,1],
type="n", col=cols[1], lwd=2,
ylim=c(0,1), main="Los Haitises",
ylab="Quasi-extinction probability", xlab="Year",
xaxt="n")
axis(1, at=seq(2020, 2120, by=10),
labels=c(2020, NA, 2040, NA, 2060, NA, 2080, NA, 2100, NA, 2120))
abline(h=seq(0, 1, by=0.1), col="gray90")
abline(v=seq(2030, 2130, by=10), col="gray90")
lines(2024:2123, pe[1, ,1], col=cols[1], lwd=2)
lines(2024:2123, pe[2, ,1], col=cols[2], lwd=2)
lines(2024:2123, pe[3, ,1], col=cols[3], lwd=2)
lines(2024:2123, pe[4, ,1], col=cols[1], lwd=2, lty=2)
lines(2024:2123, pe[5, ,1], col=cols[2], lwd=2, lty=2)
lines(2024:2123, pe[6, ,1], col=cols[3], lwd=2, lty=2)
plot(2024:2123, pe[1, ,2],
type="n", col=cols[1], lwd=2,
ylim=c(0,1), main="Punta Cana",
ylab="", xlab="Year",
xaxt="n")
axis(1, at=seq(2020, 2120, by=10),
labels=c(2020, NA, 2040, NA, 2060, NA, 2080, NA, 2100, NA, 2120))
abline(h=seq(0, 1, by=0.1), col="gray90")
abline(v=seq(2030, 2130, by=10), col="gray90")
lines(2024:2123, pe[1, ,2], col=cols[1], lwd=2)
lines(2024:2123, pe[2, ,2], col=cols[2], lwd=2)
lines(2024:2123, pe[3, ,2], col=cols[3], lwd=2)
lines(2024:2123, pe[4, ,2], col=cols[1], lwd=2, lty=2)
lines(2024:2123, pe[5, ,2], col=cols[2], lwd=2, lty=3)
lines(2024:2123, pe[6, ,2], col=cols[3], lwd=2, lty=4)
plot.new()
legend(x=-0.3, y=0.5, title="Scenarios",
legend=c("1 trans=0, nests=0", "2 trans=0, nests=10", "3 trans=0, nests=all",
"4 trans=10, nests=0", "5 trans=10, nests=10", "6 trans=10, nests=all"),
xpd=NA, col=c(cols, cols), lty=c(1,1,1,2,2,2), lwd=2,
xjust=0, yjust=0.5)
mtext("Future year", 1, outer=TRUE, line=1, adj=0.39)
mtext("Quasi-extinction probability", 2, outer=TRUE, line=1)
# dev.off()
# Time to extinction
# conditional on extinction threshold (D=X)
# D <- 3
# qextinct <- outp$Ntot <=3
# ext <- outp$extinct[,100,,]==1
# t.extinction <- apply(qextinct[])2.2 Plots of finer population segments and model diagnostics
load("C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\Ridgways IPM\\outputs\\ipm_longrun.rdata")
load("data/data.rdata")
library ('MCMCvis')
library ('coda')
library ('ggplot2')
library('reshape2')
library('bayestestR')
out <- lapply(post, as.mcmc)
# Identify chains with NAs that
# failed to initialize
NAlist <- c()
for (i in 1:length(out)){
NAlist[i] <- any (is.na(out[[i]][,1:286]) | out[[i]][,1:286]<0)
}
# Subset chains to those with good initial values
out <- out[!NAlist]
post2 <- post[!NAlist]
outp <- MCMCpstr(out, type="chains")
!NAlist## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE# default settings for plots
plt <- function(object, params,...) {
MCMCplot(object=out,
params=params,
guide_axis=TRUE,
HPD=TRUE, ci=c(80, 95), horiz=FALSE,
#ylim=c(-10,10),
...)
}Plot model estimates of demographic rates. Life Stages are abbreviated as B = breeder, NB = nonbreeder, FY = first year. First-year abundance accounts for translocated birds.
# Abundance of females at Los Haitises
par(mfrow=c(5,2))
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("NFY[",1:13, ", 1]"),
main="First-year (FY) minus hacked\n Los Haitises",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plot(2011:2023, datl$countsFY[,1],
ylab="Counts", xlab="Year", type="b")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("NF[",1:13, ", 1]"),
main="Adult nonbreeder (NB)\n Los Haitises",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plot.new()
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("NB[",1:13, ", 1]"),
main="Adult breeder (B)\n Los Haitises",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plot.new()
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("NAD[",1:13, ", 1]"),
main="Adult Breeders and Nonbreeders\n Los Haitises",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plot(2011:2023, datl$countsAdults[,1],
ylab="Counts", xlab="Year", type="b")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("Ntot[",1:13, ", 1]"),
main="All stages\n Los Haitises",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plot(2011:2023, datl$countsFY[,1]+datl$countsAdults[,1],
ylab="Counts", xlab="Year", type="b")
# Abundance of females at Punta Cana
par(mfrow=c(5,2))
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("NFY[",1:13, ", 2]"),
main="First-year (FY) plus hacked\n Punta Cana",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plot(2011:2023, datl$countsFY[,2],
ylab="Counts", xlab="Year", type="b")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("NF[",1:13, ", 2]"),
main="Adult nonbreeder (NB)\n Punta Cana",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plot.new()
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("NB[",1:13, ", 2]"),
main="Adult breeder (B)\n Punta Cana",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plot.new()
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("NAD[",1:13, ", 2]"),
main="Adult Breeders and Nonbreeders\n Punta Cana",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plot(2011:2023, datl$countsAdults[,2],
ylab="Counts", xlab="Year", type="b")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("Ntot[",1:13, ", 2]"),
main="All stages\n Punta Cana",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plot(2011:2023, datl$countsFY[,2]+datl$countsAdults[,2],
ylab="Counts", xlab="Year", type="b")
Population dynamics are determined by transitions, These transitions between stages are abbreviated as the starting life stage to the final life stage. For example a first-year recruiting to a breeder would be abbreviated as “FY to B”. I’ll list them here for convenience:
“FY to NB” is first-year to nonbreeder.
“NB to NB” is nonbreeder adult to nonbreeder adult.
“B to NB” is a breeding adult to a nonbreeder adult.
“FY to B” is first-year to breeder.
“NB to B” is nonbreeder adult to breeder adult.
“B to B” is breeder adult to breeder adult.
# Finer population segments
par(mfrow=c(4,2))
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 1, ", ", 1:13, ", 1]"),
main="Los Haitises\nFirst-years fledged",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 2, ", ", 1:13, ", 1]"),
main="Los Haitises\nFY to NB",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 3, ", ", 1:13, ", 1]"),
main="Los Haitises\nNB to NB",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 4, ", ", 1:13, ", 1]"),
main="Los Haitises\nB to NB",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 5, ", ", 1:13, ", 1]"),
main="Los Haitises\nFY to B",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 6, ", ", 1:13, ", 1]"),
main="Los Haitises\nNB to B",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 7, ", ", 1:13, ", 1]"),
main="Los Haitises\nB to B",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
par(mfrow=c(4,2))
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 1, ", ", 1:13, ", 2]"),
main="Punta Cana\nFY fledged",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 2, ", ", 1:13, ", 2]"),
main="Punta Cana\nFY to NB",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 3, ", ", 1:13, ", 2]"),
main="Punta Cana\nNB to NB",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 4, ", ", 1:13, ", 2]"),
main="Punta Cana\nB to NB",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 5, ", ", 1:13, ", 2]"),
main="Punta Cana\nFY to B",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 6, ", ", 1:13, ", 2]"),
main="Punta Cana\nNB to B",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 7, ", ", 1:13, ", 2]"),
main="Punta Cana\nB to B",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
Other parameter estimates.
# I needed to abbreviate to save plot space
# FY=first-year, NB=Nonbreeder, B=Breeder
par(mfrow=c(1,2))
plt(object=out,
params=paste0("mus[",1:8, ", 1]"),
exact=TRUE, ISB=FALSE,
ylim=c(0,1),
main="Overall means\n Los Haitises",
labels=c("FY survival", "NB survival", "B survival",
"FY to B", "NB to B", "B to NB",
"NB detection", "B detection")
)
plt(object=out,
params=paste0("mus[",1:8, ", 2]"),
exact=TRUE, ISB=FALSE,
ylim=c(0,1),
main="Overall means\n Punta Cana",
labels=c("FY survival", "NB survival", "B survival",
"FY to B", "NB to B", "B to NB",
"NB detection", "B detection"))
par(mfrow=c(1,1))
plt(object=out,
params="betas",
main= "Translocation effects",
labels=c("FY survival", "NB survival", "B survival",
"FY to B", "NB to B", "B to NB",
"NB detection", "B detection"))
# Plot effort effects
plt(object=out,
params=paste0("deltas[",1:4,"]"),
exact=TRUE, ISB=FALSE,
main= "Effort effects",
labels=c("Number of Adults\nEffort", "Number of Adults\nEffort squared",
"Number of FYs\nEffort", "Number of FYs\nEffort squared"),
ylim=c(-0.5, 1))
plt(object=out,
params=paste0("deltas[",5:8,"]"),
exact=TRUE, ISB=FALSE,
main= "Effort effects continued",
labels=c("Nonbreeder detection\nEffort", "Nonbreeder detection\nEffort sq",
"Breeder detection\nEffort", "Breeder detection\nEffort sq"),
ylim=c(-5, 5))
# Fecundity
par(mfrow=c(1,1))
plt(object=out,
params=c("lmu.prod"),
labels= c("Productivity\n(log scale)\nLos Haitises",
"Productivity\n(log scale)\nPunta Cana"))
f <- exp(outp$lmu.prod)
# Is fecundity at LHNP greater than PC
par(mfrow=c(1,1))
fdiff <- f[1,]-f[2,]
hist(fdiff, main="Productivity difference")
abline(v=0, lty=2)
# print probability of direction, similar to frequentist p-value
# so values <=0.025 and >=0.975
# Is the difference in fecundity >0 ?
mean(fdiff>0) ## [1] 1# How many times greater is fecundity
# at treated versus non-treated sites
# median(f.pred[1,]/f[1,])
# median(f.pred[2,]/f[2,])
# gamma = nest treatment effect on fecundity
par(mfrow=c(1,1))
plt(object=out,
params=c("gamma"),
main="Anti-Parasitic Fly\nTreatment Effects", ylim=c(0,3))
par(mfrow=c(1,1))
sds <- paste0("sds[", 1:9, "]")
plt(object=out, params=sds,
exact=TRUE, ISB=FALSE,
main="Temporal SDs (synchrony among sites)",
labels=c("FY survival", "NB survival", "B survival",
"FY to B", "NB to B", "B to NB",
"NB detection", "B detection",
"Productivity"))
sds2 <- paste0("sds2[", 1:9, "]")
plt(object=out, params=sds2,
exact=TRUE, ISB=FALSE,
main="Site-temporal SDs",
labels=c("FY survival", "NB survival", "B survival",
"FY to B", "NB to B", "B to NB",
"NB detection", "B detection",
"Productivity"))
# Correlations among vital rates
# Plot is messy with only a few strong correlations
ind <- 1
Rs <- R2s <- c()
for (i in 1:(nrow(outp$R)-1)){
for (j in (i+1):nrow(outp$R)){
Rs[ind] <- paste0("R[",i,", ", j, "]")
R2s[ind] <- paste0("R2[",i,", ", j, "]")
ind <- ind+1
}}
par(mfrow=c(2,1))
plt(object=out, params=Rs[1:18], exact=TRUE, ISB=FALSE,
main="Correlations btw demographic rates\n over time (synchrony)",
xlab = "Rhos", guide_lines=TRUE)
plt(object=out, params=Rs[19:36], exact=TRUE, ISB=FALSE,
main="Correlations btw demographic rates\n over time (synchrony), continued...",
xlab = "Rhos", guide_lines=TRUE)
par(mfrow=c(2,1))
plt(object=out, params=R2s[1:18], exact=TRUE, ISB=FALSE,
main="Correlations btw demographic rates\n over time and sites",
xlab = "Rhos", guide_lines=TRUE)
plt(object=out, params=R2s[19:36], exact=TRUE, ISB=FALSE,
main="Correlations btw demographic rates\n over time and sites, continued ...",
xlab = "Rhos", guide_lines=TRUE)
# Annual averages for integration into the population model
par(mfrow=c(1,1))
labs <- c(paste0("LH ",2011:2023), paste0("PC ",2011:2023))
plt(object=out, params="mn.phiFY", ylim=c(0,1),
main="First-year survival", labels = labs,
xlab = "Year", ylab= "Survival")
abline(v=13.5, lwd=2)
plt(object=out, params="mn.phiA", ylim=c(0,1),
main="Adult nonbreeder", labels = labs,
xlab = "Year", ylab= "Survival")
abline(v=13.5, lwd=2)
plt(object=out, params="mn.phiB", ylim=c(0,1),
main="Breeder", labels = labs,
xlab = "Year", ylab= "Survival")
abline(v=13.5, lwd=2)
plt(object=out, params="mn.psiFYB", ylim=c(0,1),
main="First-year to breeder", labels = labs,
xlab = "Year", ylab= "Recruitment")
abline(v=13.5, lwd=2)
plt(object=out, params="mn.psiAB", ylim=c(0,1),
main="Adult nonbreeder to breeder", labels = labs,
xlab = "Year", ylab= "Recruitment")
abline(v=13.5, lwd=2)
plt(object=out, params="mn.psiBA", ylim=c(0,1),
main="Adult breeder to nonbreeder", labels = labs,
xlab = "Year", ylab= "Recruitment")
abline(v=13.5, lwd=2)
plt(object=out, params="mn.pA", ylim=c(0,1),
main="Nonbreeder", labels = labs,
xlab = "Year", ylab= "Detection")
abline(v=13.5, lwd=2)
plt(object=out, params="mn.pB", ylim=c(0,1),
main="Breeder", labels = labs,
xlab = "Year", ylab= "Detection")
abline(v=13.5, lwd=2)
plt(object=out, params="mn.prod",
main="", labels=labs,
xlab = "Year", ylab= "Productivity")
abline(v=13.5, lwd=2)
2.3 Print parameter estimates
Parameter estimates for input into a population viability analysis.
pars1 <- c("sds", "sds2","mus", "betas",
"NFY", "NF", "NB", "Ntot",
"mn.phiFY","mn.phiA", "mn.phiB",
"mn.psiFYB", "mn.psiAB", "mn.psiBA",
"mn.pA", "mn.pB")
# Estimates for the survival model
# In this order: FY survival, NB survival, B survival,
# FY to B recruitment, NB to B recruitent, B to NB recruitment,
# Detection NB, Detection B
# Mus are means for
# mus[1, site] , where site=1 is LH and site=2 is PC
# Survival of first years = mus[1,]
# Survival of nonbreeders = mus[2,]
# Survival of breeders = mus[3,]
# Breeding propensity of first years = mus[4,]
# Breeding propensity of nonbreeders = mus[5,]
# Transition from breeder to nonbreeder = mus[6,]
# Detection probability of nonbreeders = mus[7,]
# Detection probability of breeders
# modeled as logit(probability) = lmus[x,site] + betas[x]*translocated + eps[x,t] + eta[x,s,t]
# except breeder to nonbreeder recruitment = lmus[6,site]
MCMCsummary(post2, params = "mus",
digits=3, HPD = T,
hpd_prob = 0.95, pg0= TRUE,
func=median, func_name="median")MCMCsummary(post2, params = "betas",
digits=3, HPD = T,
hpd_prob = 0.95, pg0= TRUE,
func=median, func_name="median")# Fecundity
# modeled as log(f) = lmu.f[site] + gamma*treatment + eps[x,t] + eta[x,s,t]
# log scale
MCMCsummary(post2, params = "lmu.prod",
digits=3, HPD = T,
hpd_prob = 0.95, pg0= TRUE,
func=median, func_name="median")# Effort effects on detections and resighting
MCMCsummary(post2, params = "deltas",
digits=2, HPD = T,
hpd_prob = 0.95, pg0= TRUE,
func=median, func_name="median")# Random effects for temporal synchrony and site x year.
MCMCsummary(post2, params = c(sds, sds2),
exact=TRUE, ISB=FALSE,
digits=2, HPD = T,
hpd_prob = 0.95, pg0= TRUE,
func=median, func_name="median")# Correlations among demographic rates time (synchrony)
MCMCsummary(post2, params = c(Rs, R2s),
exact=TRUE, ISB=FALSE,
digits=2, HPD = T,
hpd_prob = 0.95, pg0= TRUE,
func=median, func_name="median")# Correlations among demographic rates site x time
MCMCsummary(post2, params = R2s,
exact=TRUE, ISB=FALSE,
digits=2, HPD = T,
hpd_prob = 0.95, pg0= TRUE,
func=median, func_name="median")3.Model diagnostics
3.1 Check Goodness-of-fit
Goodness-of-fit tests provide evidence that statistical distributions adequately describe the data. Here we test fit for brood size and counts. A Bayesian p-value nearest to 0.5 suggests good fitting statistical distributions, while values near 1 or 0 suggest poor fit.
# Goodness of fit check
fit.check <- function(out, ratio=FALSE,
name.rep="f.dmape.rep",
name.obs="f.dmape.obs",
jit=100,
ind=1,
lab=""){
par(mfrow=c(1,1))
# plot mean absolute percentage error
samps <- MCMCpstr(out, "all", type="chains")
rep <- samps[name.rep][[1]][ind,]
obs <- samps[name.obs][[1]][ind,]
mx <- max(c(rep, obs))
mn <- min(c(rep, obs))
plot(jitter(obs, amount=jit),
jitter(rep, amount=jit),
main=paste0("Mean absolute percentage error\n",lab),
ylab="Discrepancy replicate values",
xlab="Discrepancy observed values",
xlim=c(mn,mx), ylim=c(mn,mx),
pch=16, cex=0.5, col="gray10")
curve(1*x, from=mn, to=mx, add=T, lty=2, lwd=2, col="blue")
bp1 <- round(mean(rep > obs),2)
loc <- ifelse(bp1 < 0.5, "topleft", "bottomright")
legend(loc, legend=bquote(p[B]==.(bp1)), bty="n", cex=2)
if (ratio==TRUE){
t.rep <- samps["tvm.rep"][[1]][ind,]
t.obs <- samps["tvm.obs"][[1]][ind,]
# plot variance/mean ratio
hist(t.rep, nclass=50,
xlab="variance/mean ", main=NA, axes=FALSE)
abline(v=t.obs, col="red")
axis(1); axis(2)
}
return(list('Bayesian p-value'=bp1))
}
# Counts of adults (nonbreeders+breeders)
fit.check(out, ratio=F,
name.rep="dmape.rep",
name.obs="dmape.obs",
ind=1,
lab="Adults(Breeder+Nonbreeder)- Poisson", jit=300)
## $`Bayesian p-value`
## [1] 0.67# Counts of first-year fledglings, ind=2
# Poisson failed fit test bp=0
# So we assigned negative binomial only for Punta Cana.
# Los Haitises excluded from neg bin because no zeroes in counts.
fit.check(out, ratio=F,
name.rep="dmape.rep",
name.obs="dmape.obs",
ind=2,
lab="First-year counts\nNeg binomial-Poisson", jit=300)
## $`Bayesian p-value`
## [1] 0.59# Productivity
fit.check(out, ratio=F,
name.rep="f.dmape.rep",
name.obs="f.dmape.obs",
ind=1,
lab="Productivity-Neg binomial", jit=300)
## $`Bayesian p-value`
## [1] 0.843.2 Examine Traceplots, Compare Posterior with Priors
Traceplots provide evidence that models adequately converged.
MCMCtrace(post2, pdf=FALSE, params= "mus", Rhat=TRUE, priors=runif(20000, 0, 1), post_zm=FALSE)
MCMCtrace(post2, pdf=FALSE, params= "betas", Rhat=TRUE, priors=runif(20000, -20, 20), post_zm=FALSE)
MCMCtrace(post2, pdf=FALSE, params= "deltas", Rhat=TRUE, priors=runif(20000, -20, 20), post_zm=FALSE)
MCMCtrace(post2, pdf=FALSE, params= "gamma", Rhat=TRUE, priors=runif(20000, -20, 20), post_zm=FALSE)
MCMCtrace(post2, pdf=FALSE, params= "sds", Rhat=TRUE, priors=rexp(20000, 1), post_zm=FALSE)
MCMCtrace(post2, pdf=FALSE, params= "sds2", Rhat=TRUE, priors=rexp(20000, 1), post_zm=FALSE)
MCMCtrace(post2, pdf=FALSE, params= "R2")
